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Introduction 


Digital  Tomosynthesis  Mammography  (DTM)  is  one  of  the  most  promising 
techniques  that  can  provide  quasi  three-dimensional  (3D)  information  of  the  breast  and 
potentially  improve  early  detection  of  breast  cancers.  The  development  of  the  DTM 
system  is  still  at  its  early  stage  and  current  image  reconstruction  methods  result  in  DTM 
images  with  strong  artifacts  due  to  the  limited-angle  data.  Moreover,  current  methods 
cannot  reconstruct  linear  attenuation  coefficients  (LACs)  accurately  due  to  scatter  and 
beam-hardening  artifacts.  The  goal  of  our  project  is  to  conduct  studies  to  develop 
advanced  DTM  reconstruction  algorithms  to  optimize  image  quality  and  minimize  image 
artifacts,  as  well  as  to  improve  quantitative  LAC  estimation  of  breast  tissues.  We  will 
focus  on  the  evaluation  of  effects  of  improved  reconstruction  algorithms  on  mass  lesions 
although  the  improved  image  quality  should  benefit  the  detection  and  diagnosis  of 
microcalcifications  as  well.  The  following  specific  aims  will  be  addressed  in  the  three- 
year  period  of  the  project:  (1)  Quantitative  comparison  of  several  algorithms  for  3D 
limited-angle  cone-beam  tomographic  problem  in  breast  tomosynthesis  mammography. 
(2)  Development  of  artifact  reduction  methods  to  suppress  multiple  sources  of  artifacts 
that  affect  image  quality  and  the  quantitative  estimation  of  LACs.  (3)  Development  of 
scatter  and  beam  hardening  correction  methods  and  further  improve  the  image  quality 
and  the  quantitative  estimation  of  LACs  in  DTM.  (4)  Evaluation  of  the  detection  accuracy 
of  mass  lesions  and  the  classification  accuracy  between  malignant  and  benign  masses 
with  and  without  the  proposed  correction  methods  using  collected  patient  data  and 
computer-aided  detection  (CADd)  and  diagnosis  (CADx)  systems  being  developed  in  our 
laboratory. 

Body 


This  is  the  first  year  annual  report  of  our  project.  In  the  project  period  (4/1/07- 
3/31/08),  we  have  conducted  a  number  of  studies  in  Digital  Tomosynthesis 
Mammography  (DTM),  including  implementation  and  comparison  of  multiple 
reconstruction  algorithms,  development  of  dedicated  breast  phantom  and  performance 
measures  for  image  quality  evaluation,  development  of  breast-shape  guided  efficient 
algorithm  for  iterative  reconstruction  method,  development  of  artifact  reduction  methods 
to  suppress  DTM  image  artifacts  from  multiple  sources,  and  investigation  of  the  impact 
of  DTM  imaging  conditions  on  image  quality. 

(A)  Implementation  and  comparison  of  multiple  algorithms  for  limited-angle  cone- 
beam  tomography  in  DTM  reconstruction,  and  development  of  dedicated  breast 
phantom  and  image  quality  evaluation  measures. 

In  DTM,  a  sequence  of  low-dose  projection-view  images  is  acquired  from  a 
limited  angular  range  while  the  total  dose  is  set  to  be  comparable  to  that  of  regular 
mammography.  The  reconstruction  of  the  3D  breast  volume  from  2D  projection  images 
therefore  represents  a  limited-angle  cone-beam  tomographic  problem.  We  have 
implemented  most  of  the  representative  algorithms  and  have  quantitatively  compared 
three  methods  based  on  dedicated  breast  phantom  studies  and  image  quality  evaluation 
measures. 

Method: 
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For  DTM  reconstruction  or  equivalently  a  limited-angle  cone-beam  tomographic 
problem,  existing  reconstruction  methods  can  be  coarsely  divided  into  four  categories: 
back-projection  algorithms,  transform  algorithms,  algebraic  reconstruction  techniques, 
and  statistical  reconstruction  algorithms.  We  have  implemented  most  of  the 
representative  algorithms,  as  listed  in  the  following: 

•  Back-projection  algorithms:  Shift-and-Add  (SAA)  method,  back-projection  (BP) 
method 

•  Transform  algorithms:  FDK  filtered  back-projection  method 

•  Algebraic  reconstruction  techniques:  algebraic  reconstruction  technique  (ART), 
simultaneous  algebraic  reconstruction  technique  (SART)  and  simultaneous  iterative 
reconstruction  technique  (SIRT) 

•  Statistical  reconstruction  algorithms:  maximum  likelihood  method  with  gradient 
algorithm  (ML-Gradient)  and  maximum  likelihood  method  with  convex  algorithm  (ML- 
Convex) 

We  have  investigated  the  impact  of  a  number  of  factors,  including  initialization, 
number  of  iteration,  relaxation  parameter  and  PV  image  accessing  strategy,  on  the 
reconstructed  image  quality  of  SART.  For  forward  projection  calculation  in  DTM,  we 
have  developed  a  ray-driven  model  in  computing  the  path-length  of  a  primary  x-ray 
intersecting  each  voxel  with  the  imaged  volume  lattice. 

We  have  developed  a  number  of  dedicated  breast  phantoms  in  image  quality 
evaluation.  The  breast  phantoms  include  simulated  masses  with  different  contrasts, 
simulated  microcalcification  clusters,  high-contrast  metal  wires  to  simulate  metal 
markers  or  biopsy  clips,  and  low-contrast  plastic  wires  to  simulate  mass  speculations. 
The  test  objects  are  distributed  in  different  regions  and  depths,  of  which  some  are 
overlapped  each  other.  Other  phantoms  include  breast-shape  slabs  consisted  of  breast- 
tissue-equivalent  materials,  i.e.  heterogeneous  mixture  of  fibroglandular-tissue- 
mimicking  material.  We  have  acquired  a  large  number  of  DTM  images  of  these 
phantoms,  as  well  as  ACR  (the  American  College  of  Radiology)  mammography 
accreditation  phantom  with  a  GE  prototype  DTM  system  installed  at  the  University  of 
Michigan  Hospital.  This  system  has  a  Csl  phosphor/a:Si  active  matrix  flat  panel  digital 
detector  with  a  pixel  size  of  0.1  mmxO.1  mm  and  the  raw  image  data  are  16  bits.  For 
tomosynthesis  imaging,  the  x-ray  tube  is  automatically  rotated  in  3°  increments  to 
acquire  projection  images  at  21  different  angles  over  a  60Q  angular  range  in  less  than  8 
seconds.  The  digital  detector  is  stationary  during  image  acquisition.  The  DTM  system 
uses  an  Rh-target/Rh-filter  x-ray  source  for  all  breast  thicknesses  and  no  anti-scatter 
grid  is  used. 

We  developed  and  implemented  a  number  of  performance  measures  in 
quantitative  evaluation  of  DTM  image  quality.  These  measures  include  the  contrast-to- 
noise  ratio,  normalized  line  profiles  of  test  objects,  depth  resolution,  intra-plane  spatial 
resolution,  inter-plane  artifact  spread  function,  noise  power  spectrum,  edge  sharpness 
measure,  and  uniformity  measure. 

Result: 

All  reconstruction  methods  that  we  implemented  can  reconstruct  the  features  in 
their  correct  layers  and  separate  superimposed  phantom  structures  along  the  Z 
direction,  provided  that  accurate  DTM  system  geometry  is  available.  Our  comparative 
study  of  multiple  representative  algorithms  suggested  that  the  two  iterative  methods, 
SART  and  ML-Convex,  both  can  provide  high-quality  reconstruction  DTM  images  in 
terms  of  the  performance  measures,  and  are  effective  in  improving  the  conspicuity  of 
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object  details  and  suppressing  inter-plane  blurring.  The  SART  method  can  provide 
comparable  tomosynthesized  image  quality  to  those  of  ML-Convex  method  but  with 
fewer  number  of  iterations. 

For  SART  method,  we  have  found  the  suitable  combinations  of  the  number  of 
iterations  and  the  relaxation  parameter  based  on  both  quantitative  analysis  of  image 
quality  measures  and  visual  comparison  of  reconstructed  image  quality  for  a  given  type 
of  cases.  For  initialization  of  the  SART  method,  both  a  constant  distribution  and  a  BP 
reconstructed  volume  have  been  evaluated.  It  is  also  found  that  the  SART  method  is 
robust  with  respect  to  PV  access  ordering.  In  addition,  for  such  a  limited-angular  range 
and  a  limited  number  of  acquired  PV  images  in  DTM,  we  have  found  that  the  order- 
subset  strategy,  which  groups  selected  part  of  the  PVs  images  together  in  updating  the 
volume  following  a  specific  order,  did  not  give  any  advantage  in  DTM  reconstruction, 
although  it  is  frequently  used  in  statistical  reconstruction  methods  for  other  full-angle 
imaging  modalities. 

Conclusion: 

The  SART  method  can  provide  high-quality  reconstruction  images  for  DTM.  The 
reconstruction  parameters  of  SART  have  to  be  further  investigated  and  optimized  when 
the  system  design  and/or  the  imaging  conditions  change. 

(B)  Development  of  efficient  algorithm  for  iterative  DTM  reconstruction  method  by 
utilization  of  breast  shape  information 

Iterative  method  such  as  SART  can  produce  good  tomosynthesized  image 
quality  compared  to  maximum  likelihood  type  algorithms.  Plowever,  the  computational 
burden  of  iterative  methods  is  the  major  concern  in  future  clinical  practice  of  DTM.  In  this 
study,  we  have  developed  methods  to  incorporate  both  2D  and  3D  breast  boundary 
information  within  the  SART  reconstruction  algorithm  to  improve  the  efficiency  of  DTM 
reconstruction.  This  will  also  address  the  problem  of  boundary  artifacts. 

Method: 

The  2D  breast  boundary  on  a  PV  image  is  segmented  by  a  breast  boundary 
detection  program  developed  in  our  laboratory  and  adapted  to  suit  DTM  application.  Two 
main  steps  were  in  the  algorithm,  including  detection  of  an  initial  breast  boundary  by 
applying  Otsu’s  thresholding  method  to  the  histogram  of  the  input  image  and 
performance  of  a  breast  boundary  tracking  procedure  based  on  the  initial  boundary  and 
gradient  information  from  both  horizontal  and  vertical  Sobel  filtering.  With  the  detected 
2D  breast  boundary  curves  on  all  PV  images,  we  developed  a  “3D  conical  trimming” 
method  to  generate  a  3D  breast  surface  within  the  imaged  volume  and  enclose  the 
compressed  breast  while  excluding  the  exterior  air  space.  This  method  restore  the 
unique  convex  hull  inside  the  3D  imaged  volume  assigned  to  the  DTM  system  by 
reconstruction  modeling  of  which  the  projections  at  different  angles  precisely  correspond 
to  the  breast  shadow  on  each  PV  image.  The  convex  hull  combined  with  the  top  and 
bottom  surfaces  of  the  breast  delineated  by  the  compression  paddle  and  the  breast 
support  plate  define  an  enclosed  breast  volume. 

We  developed  an  efficient  implementation  of  the  SART  method  based  on  the 
breast  boundary  information.  First,  when  the  imaged  volume  is  updated  by  each  PV 
image  in  the  SART  iterative  process,  the  re-projection  and  the  back-projection  of  the 
difference  between  the  calculated  data  and  the  measured  data  will  only  be  performed 
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along  those  “valid”  rays  within  the  2D  breast  area.  All  rays  outside  the  2D  breast  area 
will  be  ignored  without  any  manipulation.  Second,  after  one  complete  SART  iteration  has 
been  performed,  i.e.  all  PV  images  have  been  used  exactly  once  in  updating  the  imaged 
volume,  the  3D  breast  surface  is  used  to  mask  the  resulting  volume  and  set  the  values 
of  all  voxels  outside  the  breast  surface  to  that  of  air. 

Result: 

This  method  substantially  reduced  the  computational  effort  compared  to  full-field 
processing  by  eliminating  the  processing  of  all  rays  outside  the  breast  region  on  each 
individual  PV.  The  saving  in  computation  time  depends  on  the  breast  size  and  breast 
position  of  each  patient  case  that  will  affect  the  percentage  of  the  breast  area  with 
respect  to  the  full-field  dimension  on  PV  images.  Therefore,  the  reduction  in  computation 
time  is  more  for  small  breasts.  The  application  of  the  3D  breast  surface,  generated  from 
all  detected  2D  breast  boundaries,  to  the  reconstructed  slices  removed  all  breast 
boundary  artifacts.  With  the  2D  breast  mask  information  incorporated,  the  detector 
boundary  artifacts  outside  the  breast  have  been  eliminated.  In  addition,  the  skin-lines  in 
all  DTM  images  are  less  over-enhanced  when  the  reconstruction  is  limited  to  the  valid 
rays  within  the  2D  breast  boundary.  The  skin-thickening  artifact  of  the  original  SART 
method  is  due  to  edge  enhancement  by  the  reconstruction  algorithm:  the  inside  part  of 
breast  area  is  updated  and  the  outside  part  is  fit  to  the  background  noise  measurement, 
and  the  difference  will  be  enhanced  iteratively.  In  contrast,  with  2D  breast  mask 
information,  only  one  side  within  the  breast  boundary  has  been  updated  while  the 
outside  voxels  maintain  their  initialized  values,  therefore  the  edge  enhancement  effect  of 
the  SART  method  will  be  less  amplified  by  multiple  updates  from  the  different  PVs. 

Conclusion: 

In  this  study,  we  have  applied  the  2D  breast  boundary  information  and  the 
generated  3D  breast  surface  to  SART  reconstruction  in  breast  tomosynthesis 
mammography.  The  2D  breast  boundary  is  detected  on  the  projection  view  images  using 
a  boundary  tracking  algorithm  developed  in  our  laboratory  and  the  3D  breast  surface  is 
generated  with  a  3D  conical  trimming  method.  The  2D  breast  boundary  curves  on  PV 
images  are  used  to  restrict  the  SART  reconstruction  to  be  performed  only  inside  the 
breast  area  while  the  3D  breast  surface  is  used  to  exclude  reconstruction  artifacts 
outside  the  breast  volume.  Experimental  results  with  patient  PV  images  demonstrated 
that  the  proposed  method  can  substantially  improve  computational  efficiency  by 
eliminating  unnecessary  reconstruction  in  regions  outside  the  breast.  Both  breast 
boundary  and  detector  boundary  artifacts  can  be  effectively  removed  by  the  proposed 
method. 

(C)  Development  of  artifact  reduction  methods  to  remove  image  artifacts  from 
multiple  sources 

Due  to  the  finite  size  of  the  flat-panel  detector  and  the  limited  field  of  view,  DTM 
reconstruction  contains  artifacts  caused  by  the  truncated  projection-view  (PV)  images.  In 
this  study,  we  have  developed  methods  to  remove  two  types  of  truncation  artifacts: 
detector  boundary  discontinuity  and  underestimation  of  the  attenuation  pathlength  due  to 
the  preset  imaged  volume  in  reconstruction  modeling.  Both  of  these  truncation  artifacts 
are  apparent  mainly  at  the  image  boundary  of  DTM  slices.  They  will  obscure  the  breast 
tissue  details  near  the  boundary  of  DTMs,  potentially  affecting  the  accuracy  of  lesion 
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detection.  We  have  developed  correction  techniques  to  reduce  their  visibility  and  recover 
the  obscured  breast  tissue  details.  A  custom-built  breast  phantom  plus  a  variety  of 
image  quality  measures  were  developed  to  quantitatively  evaluate  the  improvements. 
The  proposed  methods  have  been  also  applied  to  DTM  patient  cases. 

Method: 

When  a  portion  of  the  breast  is  not  recorded  in  some  or  all  of  the  PV  images 
because  of  a  finite-size  detector,  the  truncated  projections  will  cause  artifacts  in  DTM 
reconstruction.  The  limited  field  of  view  (FOV)  results  in  unexposed  regions  in  the 
imaged  volume  and  truncation  of  the  PVs,  particularly  at  large  projection  angles.  In 
iterative  reconstruction  such  as  SART,  the  imaged  volume  is  updated  by  processing 
each  individual  PV  image,  i.e.  only  the  part  of  the  imaged  volume  within  the  cone-beam 
ray  path  of  the  current  PV  will  be  updated.  This  will  result  in  discontinuity  in  the  voxel 
values  across  the  cone-beam  path  boundary  which  will  be  enhanced  by  further  PV 
processing  in  the  same  and  subsequent  iterations.  These  artifacts,  referred  to  as 
“detector  boundary”  truncation  artifacts,  appear  as  bright  staircase-like  lines  or  bands  at 
the  two  sides  perpendicular  to  the  x-ray  source  motion  on  all  tomosynthesized  slices.  A 
local  intensity-equalization  method  was  developed  to  suppress  the  discontinuity  of  the 
reconstructed  voxel  intensity.  Specifically,  for  each  PV  image,  after  the  updating  using 
the  current  PV  is  completed,  we  replace  all  voxel  values  within  the  oblique  wedge-shape 
volume  with  some  average  values  obtained  from  their  neighborhood  region. 

A  second  source  of  truncation  artifacts  is  caused  by  the  missed  portion  of  breast 
tissue  that  is  outside  the  finite  imaged  volume  modeled  in  the  reconstruction  algorithm. 
The  missed  tissue  volume  will  cause  estimation  error  in  the  x-ray  attenuation.  Since  the 
attenuation  will  be  considered  to  have  occurred  only  within  the  shortened  pathlength,  the 
voxel  values  will  be  overestimated.  The  overestimated  attenuation  results  in  bright 
voxels  and  will  be  referred  to  as  “glaring”  artifact.  For  any  x-ray  incident  on  the  preset 
rectangular  imaged  volume  from  the  boundary  side,  it  is  assumed  that  this  ray  has 
encountered  very  large  attenuation  if  the  corresponding  detected  x-ray  intensity  is  below 
a  predetermined  threshold.  The  computed  total  attenuation  of  this  ray  will  be 
compensated  for  by  assuming  that  the  missing  volume  is  filled  with  “average”  breast 
tissue.  In  addition,  the  modified  difference  between  the  detected  data  and  the  new 
computed  data  is  normalized  by  the  compensated  total  pathlength  but  back-projected 
only  to  those  voxels  within  the  preset  imaged  volume. 

A  custom-built  breast  phantom  was  developed  for  quantitative  evaluation  of  the 
artifact  reduction.  The  phantom  contains  test  objects  including  simulated  masses  with 
different  contrast,  high-contrast  metal  wires  to  simulate  metal  markers  or  biopsy  clips, 
and  low-contrast  plastic  wires  to  simulate  spiculations.  The  test  objects  were  arranged  in 
groups  in  both  the  boundary  and  middle  regions  and  the  phantom  was  placed  at  various 
locations  with  respect  to  the  stationary  detector  such  that  the  imaged  objects  in  the 
boundary  area  were  truncated  in  varying  number  of  PVs.  Placing  the  breast  phantom  at 
the  various  locations  enables  the  investigation  of  the  effectiveness  of  the  artifact 
reduction  methods  for  different  degrees  of  PV  truncation.  In  addition  to  the  contrast-to- 
noise  ratio  (CNR)  and  the  normalized  line  profiles  of  feature  objects,  a  non-uniformity 
error  index  was  introduced  to  evaluate  the  reconstructed  image  quality  of  homogeneous 
Lucite  background  without  and  with  the  proposed  artifact  reduction  methods. 

Result: 
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Results  of  breast  phantom  study  have  demonstrated  that  the  proposed  methods 
enhanced  the  overall  reconstructed  image  quality  of  the  truncated  regions,  resulting  in 
much  smoother  reconstructed  distributions.  Substantial  increases  in  CNR  values  of  the 
test  objects  and  in  the  uniformity  of  the  reconstructed  homogeneous  Lucite  background 
were  also  obtained.  The  improved  boundary  image  areas  gained  comparable  quality  to 
that  in  the  central  region  of  view  where  no  PV  truncation  occurred.  The  artifact  reduction 
methods  do  not  affect  the  sharpness  of  the  line  profiles  of  high-density  objects.  For 
patient  DTM  reconstruction,  the  proposed  methods  effectively  suppressed  both  the 
detector  boundary  artifact  and  glaring  artifact,  which  strongly  degraded  the  visibility  of 
the  breast  structures  in  the  boundary  areas  where  the  PV  images  were  severely 
truncated.  The  local  breast  tissue  details  were  essentially  recovered. 

Conclusion: 

Both  the  detector  boundary  truncation  artifact  due  to  the  limited  detector  size  and 
the  glaring  artifact  due  to  the  limitation  of  the  imaged  volume  have  been  substantially 
reduced  and  the  overall  image  quality  of  the  truncated  boundary  regions  is  significantly 
improved.  The  recovered  breast  tissue  details  provide  useful  information  in  breast 
cancer  detection  and  diagnosis. 

(D)  Investigation  of  the  impact  of  DTM  system  and  imaging  condition  parameters 
on  the  reconstructed  image  quality 

The  development  of  DTM  system  and  the  selection  of  optimal  imaging  conditions 
are  still  in  the  early  stage.  The  image  quality  depends  on  a  lot  of  factors  including 
hardware,  system  design,  imaging  conditions,  reconstruction  algorithms  and  many  other 
parameters.  It  is  therefore  of  practical  importance  to  investigate  the  impact  of  DTM 
system  design  on  image  quality.  The  results  will  provide  insight  and  useful  information 
for  trade-offs  between  system  design  and  optimizing  image  conditions. 

Method: 

We  used  both  breast  phantom  and  developed  computer  simulation  models  to 
conduct  this  study.  The  breast  phantoms  are  similar  to  those  used  in  the 
abovementioned  studies.  For  computer  simulation,  the  same  GE  prototype  DTM  system 
was  modeled  and  forward  projection  images  of  software  phantom  were  generated.  Ultra- 
high  resolution  of  isotropic  voxel  lattice  was  used. 

We  have  investigated  a  number  of  different  imaging  conditions,  including 
different  distributions  (or  sampling)  of  PV  images  for  a  fixed  total  dose,  different 
combinations  of  number  of  PVs,  total  angular  range  and  angular  increment,  on  the 
reconstructed  image  quality.  Among  many  image  quality  measures,  we  were  particularly 
interested  in  Z-axis  resolution  (depth  resolution).  The  Z-axis  resolution  represents  the 
ability  of  the  DTM  system  to  distinguish  adjacent  objects  along  the  depth  direction  and 
therefore  is  one  of  the  most  important  advantages  provided  by  DTM  to  separate 
overlapping  breast  tissues  occurred  in  conventional  mammograms.  With  the  ultra-high 
resolution  of  the  simulation  model,  a  “point”  like  object  can  be  used  to  evaluate  Z-axis 
resolution  at  specific  location.  To  quantitatively  evaluate  the  Z-axis  resolution  for 
different  DTM  imaging  conditions,  the  vertical  line  profile  along  the  Z-axis  and  through 
the  point  object  center  was  extracted  from  the  reconstructed  volume.  The  line  profile  was 
then  normalized  to  1  at  the  maximum  value.  The  full-width-at-half-maximum  (FWHM)  of 
the  normalized  line  intensity  profile  was  used  to  evaluate  the  Z-axis  resolution.  For  a 
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given  DTM  imaging  condition,  we  also  compared  the  Z-axis  resolution  of  SART  method 
to  BP  method. 

Result: 


Preliminary  phantom  and  computer  simulation  studies  demonstrated  that:  (1)  the 
larger  the  DTM  angular  range,  the  better  the  CNR  and  ASF  of  simulated  masses;  while 
the  narrow  angular  range  gave  sharper  in-plane  edge  for  high-density  objects.  (2)  image 
sharpness  and  contrast  decreased  with  decreasing  numbers  of  PVs  used  in  the 
reconstruction.  The  interplane  artifacts  increased  with  decreasing  angular  range  of  the 
PVs.  (3)  The  Z-axis  resolution  is  relatively  independent  of  depth  within  a  typical  breast 
thickness  range  (i.e.  a  5-cm-thick  volume)  and  the  number  of  projection  views  within  a 
fixed  angular  range,  but  strongly  affected  by  tomo-angle  such  that  the  depth  resolution 
improves  with  increasing  total  tomosynthesis  angle.  The  SART  method  can  provide 
better  Z-axis  resolution  than  the  BP  method. 

Conclusion: 

DTM  system  design  and  imaging  condition  parameters  have  significant  impact  on 
reconstruction  image  quality.  Our  investigation  has  revealed  some  important  tradeoff 
between  different  parameter  combinations.  Further  work  is  underway  to  continue  the 
studies  of  these  relationships. 

Key  Research  Accomplishments 

•  Develop  breast  phantom  and  a  variety  of  figures  of  merit  for  quantitative 

evaluation  of  image  quality  and  artifact  reduction . (Task  1  (a)) 

•  Implement  multiple  three-dimensional  limited-angle  cone-beam  tomographic 

algorithms  for  DTM  reconstruction,  including  Shift-and-Add  (SAA)  method,  back- 
projection  (BP)  method,  FDK  filtered  back-projection  method,  algebraic 
reconstruction  techniques  (ART,  SART,  SIRT),  and  statistical  methods 
(maximum  likelihood  method  with  gradient  algorithm  and  maximum  likelihood 
method  with  convex  method) - (Task  1  (b)) 

•  Quantitatively  evaluate  and  compare  representative  DTM  algorithm,  including 
BP,  SART  and  ML-Convex  method;  and  investigate  the  impact  of  a  number  of 
factors,  including  initialization,  number  of  iterations,  relaxation  parameter  and  PV 

image  accessing  strategy,  on  the  reconstructed  image  quality  of  SART  . 

(Task  1  (c)  and  (d)) 

•  Develop  a  2D  breast  boundary  segmentation  method  on  PV  images  and  a  3D 

surface-generating  method  to  estimate  breast  shape  information  .  (Task  2 

(a),  (b)  and  (c)) 

•  Incorporate  both  2D  and  3D  breast  shape  information  in  iterative  DTM 

reconstruction  algorithm  to  achieve  efficient  implementation  and  reduce 
boundary  artifacts . -  (Task  2  (d)  and  (e)) 
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•  Develop  and  evaluate  artifact  reduction  methods  with  both  breast  phantom  study 
and  patient  cases  for  DTM  image  artifacts  from  multiple  sources,  including 
boundary  artifacts,  truncation  artifacts,  glaring  artifacts  and  ghost  shadow 
artifacts  of  high-density  objects - (Task  2  (e)) 

Reportable  Outcomes 

With  the  support  by  USAMRMC,  we  have  successfully  conducted  studies  in  digital 
tomosynthesis  mammography  in  this  project  year.  We  have  published  the  results  in 
peer-reviewed  journal  and  presented  the  work  in  international  conferences,  as  listed  in 
the  following. 

Peer-Reviewed  Journal  Articles: 


1.  Y.  Zhang,  H.-P.  Chan,  B.  Sahiner,  Y.-T.  Wu,  C.  Zhou,  J.  Ge,  J.  Wei,  L.  M. 
Hadjiiski,  "Application  of  boundary  detection  information  in  breast  tomosynthesis 
reconstruction",  Medical  Physics,  34(9),  3603-3613,  2007. 

Conference  Proceedings: 

1.  Y.  Zhang,  H.-P.  Chan,  Y.-T.  Wu,  B.  Sahiner,  C.  Zhou,  J.  Wei,  J.  Ge,  L.  M. 
Hadjiiski,  J.  Shi,  "Truncation  artifact  and  boundary  artifact  reduction  in  breast 
tomosynthesis  reconstruction",  Proc  SPIE  6913;  2008:  69132Y1-69132Y9. 

2.  J.  Wei,  H.-P.  Chan,  Y.  Zhang,  B.  Sahiner,  C.  Zhou,  J.  Ge,  Y.-T.  Wu,  L.  M. 
Hadjiiski,  "Classification  of  breast  masses  from  normal  tissues  on  digital 
tomosynthesis  mammography",  Proc  SPIE  6915;  2008:  6915081-6915086. 

3.  J.  Ge,  H.-P.  Chan,  B.  Sahiner,  Y.  Zhang,  J.  Wei,  L.  M.  Hadjiiski,  C.  Zhou,  Y.-T. 
Wu,  J.  Shi,  "Digital  tomosynthesis  mammography:  Improvement  of  artifact 
reduction  method  for  high-attenuation  objects  on  reconstructed  slices",  Proc 
SPIE  6913;  2008:  6913401-6913406. 

4.  H.-P.  Chan,  Y.-T.  Wu,  B.  Sahiner,  Y.  Zhang,  J.  Wei,  R.  H.  Moore,  D.  B.  Kopans, 
M.  A.  Helvie,  L.  M.  Hadjiiski,  T.  Way,  "Digital  tomosynthesis  mammography: 
comparison  of  mass  classification  using  3D  slices  and  2D  projection  views",  Proc 
SPIE  6915;  2008:  6915061-6915066.  ” 

5.  Y.  Zhang,  H.-P.  Chan,  M.  M.  Goodsitt,  A.  Schmitz,  J.  W.  Eberhard,  B.  E.  H. 
Claus,  "Investigation  of  different  PV  distributions  in  Digital  Breast  Tomosynthesis 
(DBT)  Mammography",  International  Workshop  on  Digital  Mammography 
(IWDM),  Tucson,  AZ,  July  2008,  (in  press). 


Conference  Abstracts  and  Presentations: 


1.  H.-P.  Chan,  Y.-T.  Wu,  B.  Sahiner,  Y.  Zhang,  R.  Moore,  D.  Kopans,  L.  M.  Hadjiiski 

M.  A.  Helvie,  "Digital  Breast  Tomosynthesis  Mammography:  Computerized 
Classification  of  Malignant  and  Benign  Masses",  Oral  presentation  at  the  49th 
Annual  Meeting  of  the  American  Association  of  Physicists  in  Medicine  (AAPM), 
Minneapolis,  MN.  July  22-26,  2007. 
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2. 


Y.  Zhang,  H.-P.  Chan,  B.  Sahiner,  Y.-T.  Wu,  C.  Zhou,  J.  Ge,  J.  Wei,  L.  M. 
Hadjiiski,  "Breast  tomosynthesis  reconstruction  with  simultaneous  algebraic 
reconstruction  technique  (SART):  Truncation  Artifact  and  Boundary  Artifact 
Reduction",  The  93rd  Scientific  Assembly  and  Annual  Meeting  of  RSNA, 
Chicago,  IL,  2007,  RSNA  Program  Book,  pp  316. 

3.  A.  V.  Nees,  B.  Sahiner,  Y.  Zhang,  A.  Joe,  H.-P.  Chan,  P.  Carson,  "Digital  Breast 
Tomosynthesis  vs.  Conventional  Mammography:  Early  Experience  with 
Comparison  of  Breast  Mass  Detection  and  Characterization",  The  93rd  Scientific 
Assembly  and  Annual  Meeting  of  RSNA,  Chicago,  IL,  2007,  RSNA  Program 
Book,  pp  704. 

4.  H.-P.  Chan,  Y.-T.  Wu,  B.  Sahiner,  Y.  Zhang,  R.  H.  Moore,  D.  B.  Kopans,  M.  A. 
Helvie,  L.  M.  Hadjiiski,  "Analysis  of  Mass  Characteristics  on  Digital  Breast 
Tomosynthesis  (DBT)  Mammograms:  Application  to  computer-aided  diagnosis", 
The  93rd  Scientific  Assembly  and  Annual  Meeting  of  RSNA,  Chicago,  IL,  2007, 
RSNA  Program  Book,  pp  315. 

5.  M.  A.  Helvie,  M.  A.  Roubidoux,  L.  M.  Hadjiiski,  Y.  Zhang,  P.  L.  Carson,  H.-P. 
Chan,  "Tomosynthesis  Mammography  versus  Conventional  Mammography: 
Comparison  of  Breast  Masses  Detection  and  Characterization",  The  93rd 
Scientific  Assembly  and  Annual  Meeting  of  RSNA,  Chicago,  IL,  2007,  RSNA 
Program  Book,  pp  381 . 

Conclusion 

During  the  first  year  of  the  project,  we  have  completed  a  number  of  studies  in  the 
development  of  advanced  reconstruction  algorithms  for  DTM.  We  have  implemented  and 
compared  some  representative  reconstruction  algorithms  for  limited-angle  cone-beam 
tomographic  problem  in  DTM.  The  SART  method,  a  member  of  the  family  of  algebraic 
reconstruction  techniques,  has  been  found  effective  in  providing  good  image  quality  and 
being  more  efficient  than  statistical  methods.  A  variety  of  reconstruction  parameters, 
including  initialization,  number  of  iterations,  relaxation  parameter  and  PV  accessing 
order,  were  investigated  in  an  effort  to  optimize  SART  reconstruction.  Multiple  dedicated 
breast  phantoms  and  image  quality  measures  were  developed  during  the  study  to 
facilitate  qualitative  and  quantitative  evaluations.  An  efficient  algorithm  was  developed  to 
reduce  the  computational  time  of  iterative  reconstruction  method  by  incorporating  breast 
shape  information  that  is  available  from  PV  images  and  the  DTM  geometry.  Artifact 
reduction  methods  have  been  developed  to  suppress  DTM  image  artifacts  from  multiple 
sources,  including  breast  boundary  shadows,  truncation  of  PV  images  and  the  limited 
size  of  the  modeled  imaged  volume.  The  advanced  SART  method  with  all 
abovementioned  components  can  provide  a  fast,  consistent  and  satisfactory  DTM 
reconstruction  with  improved  image  quality  and  minimized  artifacts.  In  addition,  we  have 
investigated  the  relationship  of  DTM  imaging  condition  parameters,  particularly  the 
distribution  of  PV  images,  to  reconstruction  image  quality  based  on  our  developed 
advanced  SART  method.  The  study  results  contribute  to  the  understanding  of  important 
tradeoffs  between  different  DTM  system  design  parameters. 

In  the  coming  year,  we  will  investigate  the  scatter  and  beam  hardening  effect  in 
DTM  reconstruction.  We  will  develop  artifact  reduction  methods  for  these  two  problems 
and  further  improve  DTM  reconstruction  image  quality  and  the  quantitative  estimation  of 
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x-ray  attenuation  properties  of  breast  tissue.  We  will  also  continue  our  study  in  reduction 
of  other  important  DTM  image  artifacts  that  will  benefit  the  quantitative  reconstruction  of 
DTM  imaging. 
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Digital  tomosynthesis  mammography  (DTM)  is  one  of  the  most  promising  techniques  that  can 
potentially  improve  early  detection  of  breast  cancers.  DTM  can  provide  three-dimensional  (3D) 
structural  information  by  reconstructing  the  whole  imaged  volume  from  a  sequence  of  projection- 
view  (PV)  mammograms  that  are  acquired  at  a  small  number  of  projection  angles  over  a  limited 
angular  range.  Our  previous  study  showed  that  simultaneous  algebraic  reconstruction  technique 
(SART)  can  produce  satisfactory  tomosynthesized  image  quality  compared  to  maximum  likelihood- 
type  algorithms.  To  improve  the  efficiency  of  DTM  reconstruction  and  address  the  problem  of 
boundary  artifacts,  we  have  developed  methods  to  incorporate  both  two-dimensional  (2D)  and  3D 
breast  boundary  information  within  the  SART  reconstruction  algorithm  in  this  study.  A  second 
generation  GE  prototype  tomosynthesis  mammography  system  with  a  stationary  digital  detector 
was  used  for  PV  image  acquisition  from  21  angles  in  3°  increments  over  a  ±30°  angular  range.  The 
2D  breast  boundary  curves  on  all  PV  images  were  obtained  by  automated  segmentation  and  were 
used  to  restrict  the  SART  reconstruction  to  be  performed  only  within  the  breast  volume.  The 
computation  time  of  SART  reconstruction  was  reduced  by  76.3%  and  69.9%  for  cranio-caudal  and 
mediolateral  oblique  views,  respectively,  for  the  chosen  example.  In  addition,  a  3D  conical  trim¬ 
ming  method  was  developed  in  which  the  2D  breast  boundary  curves  from  all  PVs  were  back 
projected  to  generate  the  3D  breast  surface.  This  3D  surface  was  then  used  to  eliminate  the  multiple 
breast  shadows  outside  the  breast  volume  due  to  reconstruction  by  setting  these  voxels  to  a  constant 
background  value.  Our  study  demonstrates  that,  by  using  the  2D  and  3D  breast  boundary  informa¬ 
tion,  all  breast  boundary  and  most  detector  boundary  artifacts  can  be  effectively  removed  on  all 
tomosynthesized  slices.  ©  2007  American  Association  of  Physicists  in  Medicine. 

[DOI:  10.1118/1.2761968] 
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I.  INTRODUCTION 

Currently  mammography  is  the  only  proven  method  used  for 
breast  cancer  screening.  However,  a  major  limitation  of 
mammography  is  that  the  projection  images,  recorded  on  a 
screen-film  system  (conventional  mammography)  or  by  a 
digital  detector  (digital  mammography),  only  contain  two- 
dimensional  (2D)  projection  of  three-dimensional  (3D)  ana¬ 
tomical  structures.  As  a  result,  the  sensitivity  of  breast  cancer 
detection  is  adversely  affected  by  the  camouflaging  of  can¬ 
cerous  lesions  by  dense  breast  tissues. 

Digital  tomosynthesis  mammography  (DTM)  is  one  of  the 
most  promising  techniques  that  can  potentially  improve  early 
detection  of  breast  cancers. 1-3  DTM  can  provide  3D  struc¬ 
tural  information  by  reconstructing  the  whole  imaged  vol¬ 
ume  from  a  sequence  of  projection- view  (PV)  mammograms 
that  are  acquired  at  a  small  number  of  projection  angles  over 
a  limited  angular  range.  The  total  radiation  dose  is  set  to  be 
comparable  to  that  used  in  regular  screening  mammography. 
It  has  been  demonstrated6-8  that  DTM  can  reduce  the  cam¬ 
ouflaging  effect  of  the  overlapping  hbroglandular  breast  tis¬ 
sue,  thus  improving  the  conspicuity  of  subtle  lesions.  Several 


prototype  DTM  systems  have  been  developed  based  on  digi¬ 
tal  mammography  systems5’9’10  and  preliminary  pilot  clinical 
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trials  are  being  conducted  to  evaluate  the  utility  of  DTM. 

The  reconstruction  of  the  3D  volume  of  the  compressed 
breast  from  a  sequence  of  2D  projection  images  in  DTM 
represents  a  limited-angle  cone-beam  tomographic  problem. 
Current  reconstruction  methods  for  such  a  problem  include 

•  Back-projection- type  algorithms,  including  the  shift- 
and-add  method  used  in  the  original  tomosynthesis1'14 
and  multiple  projection  method;15 

•  Fourier-transform  based  algorithms,  including  the  fil¬ 
tered  back-projection  (FBP)  method9’16'17  and  other 
transfer  function  methods,  such  as  inversion 
filtering10’18’19  and  matrix  inversion  tomosynthesis 
(MITS);20-22 

•  Algebraic  reconstruction  techniques  (ARTs),  including 
the  original  ART16,23  and  simultaneous  ART 
(SART);24’25 

•  Statistical  reconstruction  algorithms,  including  maxi¬ 
mum  likelihood  (ML)  method  with  convex 
algorithm.26’27 
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Due  to  the  limited  angular  range  employed  in  DTM  ac¬ 
quisition,  only  incomplete  projection  information  of  the  im¬ 
aged  breast  is  available  which  results  in  intra-  and  inter-  slice 
artifacts.  The  inter-slice  artifacts  come  from  the  superposi¬ 
tion  of  out-of-plane  features  on  the  in-focus  plane  which  can¬ 
not  be  removed  by  DTM  reconstruction  alone.  Many  artifact 
reduction  algorithms  have  been  investigated,  especially  those 
to  address  artifacts  of  high  attenuation  objects,  including  de- 
blurring  techniques,  order  statistics  operator,  ’  and 
voting  strategies,34  and  these  methods  have  been  incorpo¬ 
rated  into  various  reconstruction  methods. 

A  review  of  tomosynthesis  reconstruction  techniques  can 
be  found  in  Dobbins  and  Godfrey."  In  our  previous  work," 
we  have  investigated  three  representative  methods  for  DTM 
reconstruction  based  on  breast  phantom  study  and  a  second 
generation  GE  prototype  DTM  system:  the  back-projection 
method,  the  SART  and  the  maximum  likelihood  method  with 
the  convex  algorithm  (ML-Convex).  Our  comparative  study 
suggested  that  both  SART  and  ML-Convex  methods  can 
achieve  good-quality  reconstruction  and  the  SART  method 
can  provide  comparable  tomosynthesized  image  quality  to 
those  of  ML-Convex  method  but  with  fewer  number  of  it¬ 
erations. 

Two  important  aspects  should  be  considered  in  breast  to¬ 
mosynthesis  reconstruction.  Lirst,  the  demand  on  high  spatial 
resolution  results  in  a  very  high  dimensionality  of  the  inverse 
problem  and  prohibitive  computational  burden.  It  is  therefore 
important  to  investigate  methods  to  avoid  unnecessary  com- 
putation.  Some  parallel  computation  method"  and  hardware 
acceleration  techniques  have  been  recently  developed  to 
improve  the  speed  of  DTM  reconstruction.  However,  these 
methods  are  usually  dedicated  to  specific  reconstruction  al¬ 
gorithms  and  the  improvement  depends  on  the  algorithm  un¬ 
der  consideration.  Since  these  dedicated  software  or  hard¬ 
ware  cannot  be  easily  developed  for  every  reconstruction 
technique,  it  is  of  practical  importance  to  design  methods  to 
improve  the  computational  efficiency,  especially  at  the  stage 
when  the  reconstruction  algorithm  is  still  being  developed. 
Second,  the  DTM  reconstruction  is  a  severely  underdeter¬ 
mined  and  ill-posed  inverse  problem  due  to  the  limited  num¬ 
ber  of  projections  and  the  limited  angular  range  which  only 
provides  incomplete  projection  information  of  the  imaged 
breast.  As  a  result,  tomosynthesis  reconstruction  generally 
contains  strong  artifacts.  These  artifacts  are  represented  as 
repeating  ghosts  along  the  x-ray  path  directions  on  all  tomo¬ 
synthesized  slices  with  reduced  intensity.  These  artifacts  also 
appear  at  the  breast  boundary  in  which  the  pixel  intensity 
changes  abruptly  from  breast  tissue  to  air  on  PV  images.  The 
breast  boundary  artifacts,  although  easily  distinguishable 
from  breast  anatomical  features,  could  be  distracting  for  ra¬ 
diologists’  reading.  Lurthermore,  these  artifacts  could  nega¬ 
tively  affect  computerized  processing  of  DTM  images,  e.g., 
computer-aided  detection  of  mass  lesions  by  using  DTM  re- 

.  .•  ■  25,36 

construction  images. 

In  this  work,  we  developed  a  3D  conical  trimming  method 
based  on  back  projection  of  the  breast  boundaries  detected 
on  the  2D  PV  images  to  generate  the  3D  breast  surface  and 
extract  the  breast  volume.  We  investigated  the  application  of 
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Fig.  1.  Geometry  of  the  second  generation  GE  prototype  digital  tomosyn¬ 
thesis  mammography  system  used  in  this  study. 


the  2D  and  3D  breast  boundary  information  to  DTM  recon¬ 
struction  in  an  effort  to  reduce  computational  burden  and  to 
suppress  boundary  artifacts.  With  2D  breast  boundary  infor¬ 
mation  on  PV  images,  the  reconstruction  algorithm  will  only 
consider  the  x-rays  that  intercept  the  breast  volume.  This  will 
eliminate  unnecessary  computation  outside  the  breast  while 
keeping  all  useful  information  for  reconstruction.  Using  pa¬ 
tient  DTM  images,  we  demonstrate  that  the  2D  and  3D  in¬ 
formation  can  be  used  to  effectively  remove  two  boundary 
artifacts:  breast  boundary  artifacts  and  detector  boundary  ar¬ 
tifacts,  from  the  reconstructed  DTM  slices. 

II.  MATERIALS  AND  METHODS 
II. A.  Breast  tomosynthesis  system 

The  imaging  geometry  of  the  second  generation  GE  pro¬ 
totype  digital  tomosynthesis  system  for  breast  imaging  re¬ 
search  is  shown  in  Lig.  1.  The  DTM  system  has  a 
Csl  phosphor/ a :  Si  active  matrix  flat  panel  digital  detector 
with  dimensions  of  19.20  cm  X  23.04  cm  and  a  pixel  size  of 
0. 1  mm  X  0. 1  mm.  The  digital  detector  is  stationary  during 
image  acquisition.  The  DTM  system  acquires  PV  images  at 
21  different  angles  over  a  ±30°  angular  range  by  automati¬ 
cally  rotating  the  x-ray  tube  in  3°  increments.  The  distance 
from  the  x-ray  focal  spot  to  the  center  of  rotation  is  64  cm 
and  the  plane  along  which  the  x-ray  source  rotates  is  perpen¬ 
dicular  to  the  detector  surface  at  the  chest  wall  edge.  The 
focal-spot-to-detector  distance  is  66  cm.  The  image  acquisi¬ 
tion  process  takes  less  than  8  s. 

In  our  reconstruction  algorithm,  we  define  the  “imaged 
volume”  as  a  rectangular  box  having  the  same  area  as  the 
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Fig.  2.  The  use  of  vertical  and  horizontal  Sobel  filtering  in  the  breast  bound¬ 
ary  tracking  procedure. 

detector,  and  the  same  thickness  as  that  of  the  compressed 
breast  measured  by  the  DTM  system,  as  shown  in  Fig.  1. 
This  imaged  volume  is  subdivided  into  an  array  of  voxels,  of 
which  the  X  and  Y  dimensions  are  chosen  to  be  the  same  as 
the  pixel  size  of  the  detector  (0. 1  mm  X  0. 1  mm)  while  the  Z 
dimension  (the  slice  thickness)  is  chosen  to  be  1  mm  in  this 
study. 

As  described  in  our  previous  paper,  '  for  the  forward  pro¬ 
jection  model,  we  have  developed  an  algorithm  to  calculate 
the  path  length  of  a  primary  x-ray  intersecting  each  voxel 
within  the  imaged  volume  lattice,  similar  to  Siddon’s 
algorithm/  Logarithmic  transformation  is  applied  to  the  raw 
pixel  intensities  of  the  detected  image  before  reconstruction 
in  the  SART  method.  We  assume  a  monoenergetic  x-ray 
source  and  ignore  the  effects  of  scattering  and  beam  harden- 

2  38 

ing  in  this  study,  similar  to  the  approach  by  Wu  et  al.~’ 

II. B.  Detection  of  2D  breast  boundary  and  generation 
of  3D  breast  surface 

The  2D  breast  boundary  on  a  PV  image  is  segmented  by  a 
breast  boundary  detection  program  developed  in  our 
laboratory39'40  and  adapted  to  suit  DTM  application.  There 
are  two  main  steps  in  the  algorithm.  In  the  first  step,  the 
initial  breast  boundary  is  obtained  by  applying  Otsu’s  thresh¬ 
olding  method41  to  the  histogram  of  the  input  image.  The 
initial  boundary  is  only  a  rough  estimate  of  the  boundary 
location  so  that  gray  level  thresholding  is  sufficient.  In  the 
second  step,  a  breast  boundary  tracking  procedure  is  per¬ 
formed,  using  the  initial  boundary  as  a  guide,  by  extracting 
gradient  information  from  horizontal  and  vertical  Sobel  fil¬ 
tering.  The  application  of  the  Sobel  filters  is  illustrated  in 
Fig.  2.  The  right-most  pixel  of  the  initial  boundary  is  first 
determined,  as  indicated  by  point  A  in  Fig.  2.  Although  this 
point  coincides  with  the  nipple  location  in  this  figure,  it  can 
be  any  point  along  the  boundary.  The  tracking  procedure 
starts  from  point  A  and  moves  upward  and  downward,  re¬ 
spectively,  to  determine  the  final  breast  boundary.  The  verti¬ 
cal  Sobel  filter  is  used  to  enhance  the  edges  in  the  ranges 
between  A  and  B  and  between  A  and  D,  and  the  horizontal 
Sobel  filter  is  used  in  the  ranges  between  B  and  C  and  be- 
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Fig.  3.  Estimation  of  the  3D  breast  surface  from  2D  breast  boundary  curves 
by  3D  conical  trimming  of  the  volume  outside  the  breast.  Dashed  arrows 
indicate  the  potential  missed  gaps  for  which  the  breast  surface  cannot  be 
recovered. 

tween  D  and  E.  The  selection  of  either  the  vertical  or  hori¬ 
zontal  Sobel  filter  is  determined  based  on  the  slope  of  the 
tangent  to  the  current  boundary  tracking  position.  From  our 
experience,  this  relatively  simple  boundary  detection  method 
is  reliable  for  the  DTM  PV  images  we  have  processed  so  far. 
However,  the  focus  of  this  work  is  not  on  the  development  of 
image  segmentation  methods.  Our  breast  surface  reconstruc¬ 
tion  algorithm  described  below  does  not  depend  on  the 
boundary  detection  method  so  that  any  satisfactory  breast 
boundary  curves  for  DTM  PV  images  obtained  with  our 
method  or  other  segmentation  techniques  can  be  used. 

With  the  detected  2D  breast  boundary  curves  on  all  PV 
images,  we  can  generate  a  3D  breast  surface  within  the  im¬ 
aged  volume  to  enclose  the  compressed  breast  while  exclud¬ 
ing  the  exterior  air  space.  In  principle,  there  exists  a  unique 
convex  hull  inside  the  3D  imaged  volume  for  the  DTM  sys¬ 
tem  of  which  the  projections  at  different  angles  precisely 
correspond  to  the  breast  shadow  on  each  PV  image.  To  re¬ 
cover  this  convex  hull,  we  implement  a  “3D  conical  trim¬ 
ming”  process.  This  process  is  schematically  depicted  in  Fig. 
3  where  only  three  x-ray  source  positions  are  drawn  for  il¬ 
lustration.  Specifically,  for  a  given  PV  image,  the  breast 
boundary  curve  is  back  projected  into  the  imaged  volume 
and  all  “air  volume”  outside  the  generalized  cone  surface 
formed  by  the  x-ray  source  and  the  breast  boundary  is 
trimmed  off.  This  process  will  repeat  for  all  PV  images.  As 
the  x-ray  source  moves  from  one  to  the  next  position,  the 
back-projected  breast  surface  will  trim  off  more  air-volume 
in  addition  to  the  preceding  ones.  After  all  PV  breast  bound¬ 
ary  curves  have  been  used  exactly  once,  the  surface  of  the 
remaining  volume  becomes  the  desired  convex  hull  and  its 
projection  with  respect  to  each  x-ray  source  position  corre¬ 
sponds  to  exactly  the  breast  shadow.  The  convex  hull  com¬ 
bined  with  the  top  and  bottom  surfaces  of  the  breast  delin¬ 
eated  by  the  compression  paddle  and  the  breast  support  plate 
define  an  enclosed  breast  volume. 

Note  that  for  any  x-ray  source  position,  some  parts  of  the 
imaged  volume  will  not  be  included  in  the  cone-beam  ray 
path,  and  at  larger  projection  angles,  some  parts  of  the  pro¬ 
jected  image  will  be  cut  off  due  to  the  limited  detector  size. 
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Fig.  4.  PV-count  maps  of  three  selected  depths:  (a)  5.0  cm,  (b)  2.5  cm,  and 
(c)  1.0  cm  above  the  breast  support  plate.  The  darker  pixel  intensity  repre¬ 
sents  a  fewer  number  of  effective  PVs.  For  all  three  slices,  the  top-most  area 
contains  only  eight  PVs  and  the  bottom-most  area  contains  only  seven  PVs. 
This  number  gradually  increases  to  2 1  from  both  ends  to  the  middle  region. 


This  problem  is  accentuated  in  a  DTM  system  with  a  station¬ 
ary  detector.  For  different  x-ray  source  positions,  the  unex¬ 
posed  parts  of  the  imaged  volume  change.  The  overall  effect 
is  that  some  middle  part  of  the  imaged  volume  is  completely 
inside  the  x-ray  path  at  all  times,  i.e.,  exposed  by  all  21  PVs; 
and  for  other  parts,  it  is  exposed  by  only  a  subset  of  the  21 
PVs.  Figure  4  shows  examples  of  the  PV-count  maps,  in 
which  the  intensity  of  a  given  pixel  is  proportional  to  the 
number  of  PVs  contributing  radiation  to  this  pixel.  The  PV- 
count  maps  of  three  tomosynthesized  slices  at  depths,  5.0, 
2.5,  and  1.0  cm,  above  the  breast  support  plate,  are  shown.  It 
can  be  seen  that  all  PVs  contribute  to  the  middle  region  of 
each  tomosynthesized  slice,  and  this  region  grows  larger  as 
the  distance  from  the  support  plate  decreases.  The  number  of 
steps  on  each  side  is  determined  by  the  imaging  geometry  of 
the  DTM  system,  and  is  equal  to  the  number  of  PVs  of  which 
the  x-ray  field  boundary  intersects  the  imaged  volume  on  that 
side.  Furthermore,  for  any  tomosynthesized  slice,  the  closer 
the  region  is  to  the  imaged  volume  boundary,  the  fewer  the 
number  of  PV s  will  expose  the  region,  as  apparent  from  the 
DTM  geometry.  In  this  specific  example  that  we  acquired 
with  our  GE  prototype  DTM  system,  only  eight  PVs  can 
expose  the  uppermost  while  seven  for  the  lowest  boundary 
area  on  all  tomosynthesized  slices.  Therefore  13  and  14 
steps,  respectively,  can  be  seen  on  the  upper  and  lower  sides 
of  the  PV-count  map.  The  unequal  PV-count  numbers  on  the 
upper  and  lower  boundary  areas  may  be  attributed  to  the 
slightly  asymmetric  angular  positions  of  the  x-ray  source 
relative  to  the  central  location. 

According  to  the  3D  conical  trimming  process,  for  the 
fully  exposed  parts  of  the  imaged  volume  (in  terms  of  PVs), 
the  3D  breast  surface  tangential  to  the  x  rays  can  be  gener¬ 
ated  by  the  intersection  of  all  back-projected  2D  breast 
boundary  curves.  Thus,  a  voxel  can  be  considered  to  be  in¬ 
side  the  breast  volume  if  and  only  if  its  projection  is  within 
the  breast  boundary  in  all  PVs.  In  contrast,  in  any  partially 
exposed  parts  of  the  imaged  volume,  a  voxel  is  considered  to 
be  inside  the  breast  volume  as  long  as  its  projection  is  within 
the  breast  boundary  of  any  PV.  This  process  ensures  that  any 


potential  breast  volume  which  may  be  contained  only  in 
some  of  the  PV  images  due  to  the  detector  cutoff  will  be 
included.  All  voxels  in  the  air  volume  above  and  below  the 
breast  volume  will  be  excluded  by  the  breast  boundaries  de¬ 
fined  by  the  compression  paddle  and  the  breast  support  plate. 


II.C.  Simultaneous  algebraic  reconstruction  technique 
(SART)  and  boundary  artifacts 

SART  method-4'42  is  an  iterative  reconstruction  technique. 
It  was  proposed  as  a  refinement  of  the  original  ART.16  4 '44  In 
SART  method,  an  update  of  the  3D  distribution  of  the  attenu¬ 
ation  properties  of  the  imaged  volume  is  performed  based  on 
each  PV.  The  calculated  current  projection  will  be  compared 
to  the  actual  detection  data  and  the  difference  will  be  back 
projected  and  added  to  the  imaged  volume  distribution.  The 
updating  is  performed  after  all  rays  in  one  PV  have  been 
processed,  resulting  in  a  “block-iterative”  strategy.  In  con¬ 
trast,  ART  method  updates  the  distribution  based  on  a  ray- 
by-ray  manner.  The  number  of  updates  in  one  complete  it¬ 
eration  of  SART  is  equal  to  the  number  of  PVs. 

Let  the  whole  image  volume  be  subdivided  into  J  voxels, 
and  the  linear  attenuation  coefficient  for  the  jth  voxel  be 
denoted  by  Xj,  1  <  7;  the  digital  area  detector  contains  I 

elements,  and  the  /th  ray,  I  <  z  <  /,  is  defined  as  a  line  seg¬ 
ment  starting  from  the  point  x-ray  source  location  to  the 
center  of  the  /th  detector  element.  The  number  of  rays  is 
equal  to  the  number  of  detector  elements  assuming  one  ray  is 
traced  for  each  element.  The  path  length  of  the  /th  ray  going 
through  the  /th  voxel  in  the  nth  x-ray  tube  location  (projec¬ 
tion  view)  is  denoted  by  ajj  n,  resulting  in  a  matrix-vector 
form  of  the  projection  model  as 

A„x  =  y„,  (1) 


where  A„  is  the  projection  matrix  for  the  nth  PV  with  uVjn  as 
its  (Z  ,y)th  element;  y„  is  the  corresponding  vector  of  the  pro¬ 
jection  data,  which  can  be  derived  from  the  pixel  values  of 
the  detected  image,  1  <  n  <  N,  where  N  is  the  total  number  of 
PVs.  The  /th  projection  value,  yin,  is  proportional  to  the 
logarithmic  transform  of  the  ratio  of  the  incident  intensity 
(70  „)  and  the  transmitted  intensity  (/,■„)  of  the  /th  ray 


yi,n  =  k  In  : 


(2) 


We  stack  the  projection  model  in  Eq.  (1)  for  all  PVs  to¬ 
gether  as 


A-l 

yt 

X  = 

_Aa,_ 

_y  n. 

(3) 


With  the  linear  system  model  in  Eq.  (3),  the  SART 
method  can  be  formulated  in  a  matrix-vector  form  as 


x"  =  x"-1  +  X  •  M„AX(y„  -  A„x"-‘) ,  (4) 

where  M„  e  RJXJ  and  W„  e  RIXI  are  diagonal  matrices  con¬ 
taining  diagonal  elements  as  [M,,]^ = (S[= ,  a,;n)- 1  and 
[ W„] a = (24= ,  aij  n) ~ 1 ,  respectively.  One  can  see  that  M„  and 
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W„  contain  inverse  “column  sum”  and  “row  sum”  of  the 
projection  matrix  A„  that  correspond  to  the  total  path  length 
of  all  x-rays  through  one  imaged  voxel  and  the  total  path 
length  of  one  x-ray  through  the  entire  imaged  volume,  re¬ 
spectively.  The  updating  strategy  of  SART  is  easy  to  under¬ 
stand  based  on  Eq.  (4)  as  follows:  the  calculated  projection 
A„x"_I  is  compared  to  the  current  PV  data  y„;  the  difference 
(y„-A„x"_1)  is  normalized  by  the  total  path  length  of  each 
x-ray  W„,  yielding  the  estimated  weighted  average  differ¬ 
ence  in  the  linear  attenuation  coefficient  of  the  voxels  along 
the  ray.  This  average  difference  is  back  projected  to  the  im¬ 
aged  volume  by  Ajt,  and  then  for  any  given  voxel,  the  dif¬ 
ferences  are  further  normalized  by  the  total  path  lengths  of 
all  rays  traveling  through  this  voxel  by  M„.  The  voxel  is 
updated  as  the  sum  of  the  current  estimation  x""1  and  the 
weighted  average  difference  scaled  by  a  relaxation  parameter 
X.  The  updating  operation  starts  from  an  initial  distribution 
x°  and  is  repeated  for  all  PV  images  using  the  corresponding 
projection  matrices.  One  complete  iteration  is  defined  as  the 
process  that  all  PV  images  have  been  sequentially  used  ex- 
actly  once.  In  our  previous  work,  we  experimentally  dem¬ 
onstrated  that  SART  method  can  provide  satisfactory  tomo- 
synthesized  image  quality  with  one  complete  iteration  using 
all  PV  images  and  \  =  0.5. 

In  this  work,  we  focus  on  two  types  of  boundary  artifacts: 
one  from  the  breast  boundary  and  the  other  from  the  detector 
boundary.  Breast  boundary  artifacts  are  caused  by  the  large 
difference  in  the  pixel  intensities  in  the  two  areas  inside  and 
outside  the  2D  breast  boundary,  corresponding  to  projection 
of  overlapping  breast  tissues  and  air  volume,  respectively. 
Due  to  the  uncertainty  in  the  true  3D  location  of  any  feature 
within  the  imaged  volume  based  on  a  limited  number  of  PV 
images  over  a  limited  acquisition  angular  range,  when  back 
projected  by  DTM  reconstruction  methods,  these  boundaries 
in  the  PVs  will  result  in  repeating  ghost  artifacts  of  similar 
breast  shape  but  reduced  intensity  on  all  tomosynthesized 
image  slices. 

Detector  boundary  artifacts  are  caused  by  the  DTM  sys¬ 
tem  geometry  in  which  a  stationary  finite-size  detector  ac¬ 
quires  PV  images  of  the  breast  from  moving  x-ray  source 
positions,  as  discussed  in  Sec.  II  B.  In  every  PV  position,  a 
part  of  the  imaged  volume  is  not  exposed  by  the  cone-beam 
ray  path  and  recorded  by  the  detector,  as  shown  in  Fig.  5(a). 
The  unexposed  part  of  the  imaged  volume  varies  as  the  x-ray 
source  moves.  When  we  reconstruct  the  imaged  volume  us- 
ing  the  SART  method,"  the  voxel  values  will  be  updated 
iteratively  by  processing  each  individual  PV  image.  For  a 
given  PV,  only  the  part  of  the  imaged  volume  within  the 
current  field  of  view  will  be  updated.  Since  the  voxels  within 
the  PV  image  boundary  (i.e.,  the  detector  boundary)  will  be 
changed  whereas  the  neighboring  voxels  outside  will  main¬ 
tain  their  previous  values,  it  will  cause  a  discontinuity  in  the 
voxel  values  across  the  image  boundary.  Similar  situation 
occurs  when  the  breast  image  exceeds  the  detector  area  and 
the  real  breast  boundary  is  outside  the  detector  boundary,  as 
schematically  shown  in  Fig.  5(b).  This  situation  occurs  most 
frequently  on  the  PV  images  associated  with  large  projection 


X-ray  source 


(a) 


X-ray  source 


(b) 

Fig.  5.  Breast  boundary  artifacts  (indicated  by  open  stars)  and  detector 
boundary  artifacts  (indicated  by  open  crosses)  appear  on  all  tomosythesized 
slices  along  the  x-ray  path  for  one  x-ray  source  position.  The  imaged  vol¬ 
ume  is  defined  as  the  volume  enclosed  by  the  rectangular  box  while  the 
breast  volume  is  that  enclosed  by  the  ellipse.  This  example  shows  the  situ¬ 
ations  (a)  that  the  breast  boundary  is  completely  contained  within  the  detec¬ 
tor  boundary;  (b)  that  one  side  of  the  breast  boundary  is  cut  off  by  the 
detector  boundary  due  to  both  the  large  breast  size  and  the  wide  projection 
angle. 


angles,  in  the  cases  of  very  large  breast  sizes,  or  at  the  top 
part  of  the  pectoral  muscle  on  most  MFO  views. 

As  an  example,  Fig.  6(a)  shows  a  patient  CC-view  DTM 
slice  reconstructed  by  the  SART  method  with  one  iteration. 
The  breast  boundary  artifacts  were  clearly  seen  as  repeated 
breast-shaped  shadows  and  the  detector  boundary  artifacts  as 
multiple  horizontal  lines  in  both  the  top  and  the  bottom  ar¬ 
eas. 

The  breast  boundary  artifacts  generally  will  be  enhanced, 
rather  than  reduced,  by  the  SART  method  as  the  number  of 
iteration  increases  until  a  convergent  solution  is  reached. 
This  is  because  the  SART  method,  similar  to  other  iterative 
reconstruction  algorithms,  acts  like  an  unsharp-masking 
(bandpass)  filter  in  updating  the  imaged  volume  and  thus 
favors  edge  enhancement.  The  convergent  solution,  however, 
will  generally  be  very  noisy  due  to  the  severe  ill-posed  na¬ 
ture  of  the  inverse  problem. 

For  detector  boundary  artifacts,  they  may  be  enhanced  or 
reduced  in  the  process  of  approaching  the  convergent  solu- 
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Fig.  6.  One  selected  tomosynthesized  slice  obtained  by  SART  reconstruc¬ 
tion  method  with  (a)  one,  (b)  two,  and  (c)  three  iterations,  without  using  any 
breast  boundary  information.  The  breast  boundary  artifacts  were  clearly  ob¬ 
served  as  repeated  breast-shaped  shadows  while  the  detector  boundary  arti¬ 
facts  as  multiple  horizontal  lines  in  both  the  top  and  the  bottom  areas.  The 
same  window  level  and  window  width  were  used  for  display  of  the  images. 
The  artifacts  did  not  diminish  as  the  number  of  iterations  increased.  Open 
rectangles  in  (a)  indicate  the  regions  of  interest  for  which  the  line  profiles 
are  compared  in  Fig.  7. 


tion,  depending  on  the  access  order  of  the  PVs  in  the  itera¬ 
tive  reconstruction  process.  To  demonstrate  these  effects,  we 
applied  the  SART  reconstruction  method  to  the  patient  CC- 
view  DTM  up  to  three  iterations,  and  presented  the  result  of 
the  same  tomosynthesized  slice  in  Figs.  6(b)  and  6(c),  re¬ 
spectively.  A  constant  relaxation  parameter  of  0.5  was  used 
for  all  three  SART  iterations.  Four  regions  of  interest,  de¬ 
picted  by  open  rectangles  in  Fig.  6(a),  were  selected  for  com¬ 
parison.  Regions  A  and  D  are  located  in  the  detector  bound¬ 
ary  artifact  area  while  regions  B  and  C  in  the  breast 
boundary  artifact  area.  The  line  profiles  in  these  four  regions 
are  shown  in  Fig.  7.  The  enhancement  of  both  types  of 
boundary  artifacts  by  SART  iterations  can  be  observed 
clearly  in  regions  A,  B  and  C.  In  region  D,  the  detector 
boundary  artifacts  were  more  prominent  than  those  in  region 
A  at  the  upper  part.  This  difference  can  be  attributed  to  the 
access  order  of  the  PV  images  employed  in  SART  recon¬ 
struction  of  this  example,  sequentially  from  that  with  the 
x-ray  source  at  the  top  side  to  the  bottom  side.45  The  upper 
part  of  the  imaged  volume  were  covered  completely  and  thus 
updated  uniformly  by  the  first  eight  PVs  at  the  beginning  of 
the  iteration.  The  voxels  of  the  entire  region  therefore  at¬ 
tained  reasonable  values  before  the  subsequent  PVs  that  par¬ 
tially  covered  the  upper  part  of  the  imaged  volume  started  to 
generate  the  detector  boundary  artifacts.  On  the  contrary,  the 
artifacts  at  the  lower  part  were  generated  starting  from  the 
first  PV.  The  large  difference  between  the  initialization  con¬ 
stant  values  outside  the  PV  image  and  the  updated  voxel 
values  within  the  PV  image  resulted  in  large  discontinuity  in 
the  voxel  values  across  the  detector  boundary.  Although  the 


(d)  Pixel  Index  (0.1  mm/pixel) 


Fig.  7.  comparison  of  line  profiles  extracted  from  selected  detector  bound¬ 
ary  artifact  areas  (regions  A  and  D)  and  breast  boundary  artifact  areas  (re¬ 
gions  B  and  C)  of  SART  results  with  one  (solid  line),  two  (dotted  line)  and 
three  iterations  (dashed  line).  The  locations  of  these  selected  regions  are 
indicated  by  open  rectangle  in  Fig.  6(a).  Pixel  index  always  starts  from  the 
end  that  is  close  to  the  bottom  part  of  the  tomosynthesized  slice. 
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subsequent  iterations  still  enhanced  the  edges,  their  dominant 
effect  was  to  bring  the  average  voxel  values  to  the  level  of 
those  of  air  attenuation.  The  overall  result  was  the  apparent 
reduction  of  edge  artifacts  with  further  iterations  in  region  D, 
opposite  to  those  observed  in  the  other  regions. 

To  verify  these,  we  continued  the  SART  reconstruction  up 
to  ten  iterations  and  observed  that  the  iteration  appeared  to 
have  reached  a  convergent  solution  for  this  DTM  case  (re¬ 
sults  not  shown  here).  The  reconstructed  background  was 
close  to  the  local  background  level  and  the  boundary  artifacts 
became  more  symmetric  in  the  upper  and  lower  parts  of  the 
slices.  The  line  profiles  of  regions  A  and  D  in  the  convergent 
result  exhibited  similar  detector  boundary  artifacts  which  did 
not  vanish  with  further  iterations.  Note  that,  because  the 
large  number  of  iterations  causes  over  enhancement  of  the 
breast  structures  and  excessive  computation  time,  the  recon¬ 
struction  will  generally  terminate  long  before  it  reaches  the 

25 

convergent  solution.  We  have  found  in  our  previous  study" 
that  good-quality  patient  DTMs  can  be  achieved  after  only 
one  iteration  with  SART. 


II. D.  Application  of  the  2D  and  3D  breast  boundary 
information  to  SART  reconstruction 

We  evaluated  the  application  of  breast  boundary  informa¬ 
tion  to  DTM  reconstruction  in  two  ways  with  the  SART 
method.  First,  when  we  update  the  imaged  volume  in  the 
iterative  process,  the  re-projection  and  the  back-projection  of 
the  difference  between  the  calculated  data  and  the  measured 
data  will  only  be  performed  along  those  “valid”  rays  within 
the  2D  breast  area.  We  will  not  use  any  other  rays  and  thus 
the  voxels  in  the  imaged  volume  associated  with  these  “in¬ 
valid”  rays  are  not  updated  and  stay  at  the  current  values. 
Second,  after  one  complete  SART  iteration  has  been  per¬ 
formed,  we  use  the  3D  breast  surface  to  mask  the  resulting 
volume  and  set  the  values  of  all  voxels  outside  the  breast 
surface  to  that  of  air. 

In  terms  of  Eq.  (4),  the  use  of  the  2D  breast  curves  and 
the  3D  breast  surface  can  each  be  expressed  as  a  linear  op¬ 
erator  in  matrix  form,  denoted  by  diagonal  matrices  P„ 
e  R,XI  and  Q  e  RJXJ,  respectively,  where  the  diagonal  ele¬ 
ment  of  P„  is  1  if  the  corresponding  detector  element  is 
marked  as  breast  area  in  the  nth  PV,  and  0  as  air  area;  simi¬ 
larly,  the  diagonal  element  of  Q  is  1  if  the  corresponding 
voxel  is  marked  as  breast  volume,  and  0  as  air  volume.  By 
using  these  linear  operators,  the  SART  method  can  be  modi¬ 
fied  as  follows: 

x"  =  x""1  +  X  •  M„A^W„P„(y„  - 

(5) 

xN=QxN. 

Note  the  operator  Q  is  applied  after  one  SART  iteration  is 
completed. 

Our  proposed  method  reduces  the  computational  effort  by 
eliminating  the  processing  of  all  rays  outside  the  breast  re¬ 
gion  on  each  individual  PV.  Thus  the  improved  efficiency 
ratio  can  be  defined  as  the  summation  of  the  background 


(a)  (b)  (c) 


Fig.  8.  Projection  view  images  of  a  DTM  in  MLO  view:  (a)  PV  at  -30°,  (b) 
PV  at  0°,  and  (c)  PV  at  +30°.  The  x-ray  source  moved  in  the  vertical 
direction  relative  to  the  images  shown.  The  detected  breast  boundary  was 
indicated  by  the  white  solid  line  on  each  PV  image 

areas  outside  the  detected  breast  boundary  on  all  PV  images 
divided  by  the  product  of  the  detector  area  and  the  number  of 
PVs. 

III.  RESULTS 

The  proposed  method  for  artifact  reduction  has  been 
evaluated  by  using  patient  PV  images  acquired  with  the  sec¬ 
ond  generation  GE  prototype  DTM  system.  Institutional  Re¬ 
view  Board  approval  and  informed  consent  were  obtained 
for  collection  of  the  patient  cases.  DTM  images  were  ac¬ 
quired  in  two  views:  a  cranio-caudal  (CC)  view  and  a  me- 
diolateral  oblique  (MLO)  view.  In  the  rest  of  the  paper,  only 
image  results  of  MLO  view  were  shown  without  loss  of  gen¬ 
erality.  CC  view  results  were  similar  to  those  of  MLO  views 
on  the  side  without  the  pectoral  muscle.  The  No.  1,  No.  11 
and  No.  21  raw  PV  images  of  MLO  view,  corresponding  to 
the  projection  angles  -30°,  0°,  +30°,  are  shown  in  Fig.  8. 
The  detected  breast  boundaries  were  indicated  by  solid  white 
lines.  The  x-ray  source  moved  in  the  vertical  direction  rela¬ 
tive  to  the  images  shown.  Figure  9  shows  the  generated  3D 


Fig.  9.  The  3D  breast  surface  generated  from  2D  breast  boundary  curves  for 
MLO  view  by  using  the  3D  conical  trimming  method. 
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Fig.  10.  Three  selected  tomosynthesized  slice  obtained  by  SART  reconstruction  for  MLO  view:  (a),  (d),  and  (g)  3.9  cm,  (b),  (e)  and  (h)  3.3  cm,  (c),  (f),  and 
(i)  2.3  cm  above  the  breast  support  plate;  (a)-(c)  reconstructed  slices  without  using  any  breast  boundary  information,  (d) — (f)  reconstructed  slices  with  using 
only  2D  breast  mask  information  (g)-(i)  reconstructed  slices  with  both  the  2D  and  3D  breast  boundary  information.  Note  that  with  2D  breast  mask,  the 
detector  boundary  artifacts  at  the  bottom  have  been  eliminated  while  the  breast  boundary  artifacts  have  been  totally  removed  after  application  of  the  3D  breast 
shape.  The  skin  lines  are  also  less  over-enhanced.  The  detector  boundary  artifacts  at  the  top  of  the  slice  do  not  change  because  the  tissue  extended  beyond  the 
detector  boundary.  SART  results  with  different  methods  ((a)— (c),  (d)-(f),  and  (g)— (i))  are  displayed  with  different  window  width  and  level  to  achieve  visually 
comparable  contrast  of  the  boundary  artifacts. 


breast  surface  from  the  MLO  view  by  using  the  proposed  3D 
conical  trimming  method. 

Figure  10  shows  examples  of  MLO-view  DTM  slices  re¬ 
constructed  with  the  SART  method  using  one  iteration  and  a 
relaxation  parameter  of  0.5.  A  constant  distribution  of  voxel 
values  was  used  as  initialization.  The  DTM  images  recon¬ 
structed  without  using  boundary  information,  as  shown  in 
(a)-(c);  with  only  2D  breast  mask  information,  as  shown  in 
(d)-(f);  and  with  both  2D  mask  and  3D  breast  shape  infor¬ 
mation,  as  shown  in  (g)-(i),  were  compared.  The  results  of 
SART  method  with  only  3D  breast  shape  information  should 
be  the  same  as  those  with  both  2D  and  3D  breast  informa¬ 
tion,  but  one  will  not  gain  computational  efficiency.  The 
three  chosen  slices  are  3.9,  3.3,  and  2.3  cm  above  the  breast 
support  plate  and  the  thickness  of  the  imaged  breast  was 


measured  by  the  DTM  system  as  5.4  cm.  While  the  choice  of 
these  tomosynthesized  slices  was  somewhat  arbitrary,  we  at¬ 
tempt  to  choose  slices  that  can  clearly  demonstrate  the  dif¬ 
ference  between  the  results  with  and  without  using  breast 
boundary  information.  Note  that  in  Figs.  10(d)-10(f),  the 
background  voxels  retained  the  initialization  constant  as  their 
values,  which  was  chosen  to  be  the  average  linear  attenuation 
coefficient  of  normal  breast  tissues  and  was  higher  than  those 
of  the  fatty  tissue.  In  contrast,  the  exterior  region  outside  the 
breast  in  Figs.  10(g)— 10(i)  were  set  to  have  values  of  zero, 
which  was  approximately  the  attenuation  coefficient  of  air, 
after  application  of  the  3D  breast  shape  information. 

Without  using  any  breast  boundary  information,  the  detec¬ 
tor  boundary  artifacts  were  clearly  observed  as  multiple  hori¬ 
zontal  lines  especially  at  the  bottom  of  all  three  tomosynthe- 
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sized  slices.  At  the  top  part  of  these  slices,  there  were  also 
detector  boundary  artifacts,  as  indicated  by  the  PV-count 
map  shown  in  Fig.  4,  but  the  artifacts  were  much  weaker 
than  those  at  the  bottom  as  a  result  of  the  PV  access  order, 
discussed  above.  The  breast  boundary  artifacts,  represented 
as  repeated  breast-shaped  ghosts  around  the  boundary,  were 
clearly  observed  around  the  true  breast  boundary  in-focus  on 
each  slice.  The  number  of  ghost  shadows  is  20  if  the  breast 
boundary  is  imaged  in  all  21  PVs  because  one  boundary  is 
in-focus  and  appears  as  real  breast  boundary  at  that  depth. 

The  application  of  the  generated  3D  breast  surface  to  the 
reconstructed  slices  removed  all  breast  boundary  artifacts. 
With  the  2D  breast  mask  information  incorporated,  the  de¬ 
tector  boundary  artifacts  outside  the  breast  have  been  elimi¬ 
nated.  More  importantly,  the  skin  lines  in  all  DTM  images 
are  less  over  enhanced  when  the  reconstruction  is  limited  to 
the  valid  rays  within  the  2D  breast  boundary.  The  over¬ 
enhanced  skin  lines  may  create  skin-thickening  artifacts.  The 
remaining  detector  artifacts  within  the  generated  breast  sur¬ 
face,  as  shown  at  the  top  part  of  the  slices,  are  due  to  the  fact 
that  the  pectoral  muscle  was  projected  beyond  the  detector 
boundary  at  the  large  tomosynthesis  angles,  as  shown  in  Fig. 
8. 

In  this  case,  the  computation  time  was  reduced  by  76.3% 
and  69.9%  for  the  CC  and  MLO  views,  respectively,  by  us¬ 
ing  the  2D  breast  boundary  masks  for  all  PVs. 

IV.  DISCUSSION 

We  have  investigated  the  application  of  the  2D  and  3D 
breast  boundary  information  to  DTM  reconstruction  with  the 
SART  method.  Preliminary  results  show  that  this  method  can 
substantially  reduce  the  computation  time  and  effectively  re¬ 
move  the  breast  and  detector  boundary  artifacts  on  tomosyn- 
thesized  slices  by  employing  only  projection  information 
from  the  rays  intersecting  the  breast  volume  for  reconstruc¬ 
tion.  When  the  breast  boundary  is  completely  inside  the  de¬ 
tector  boundary,  the  use  of  valid  rays  only  within  the  breast 
areas  for  reconstruction  excludes  the  transition  regions  in  the 
imaged  volume  that  are  covered  by  different  number  of  PVs 
and  enhanced  during  updating  in  subsequent  PV  processing. 
The  detector  boundary  artifacts  can  thus  be  totally  removed. 
When  the  projected  breast  image  exceeds  the  detector  area, 
the  detector  boundary  artifacts  outside  the  breast  volume  can 
be  removed  but  those  inside  remain  where  the  breast  PV 
images  are  cut  off  by  the  detector  boundary.  After  the  SART 
reconstruction,  breast  boundary  artifacts  were  removed  by 
trimming  off  the  air  volume  outside  the  3D  breast  surface. 
The  investigation  of  methods  for  reduction  of  the  remaining 
artifacts  will  be  undertaken  in  the  future. 

The  2D  and  3D  breast  boundary  information  can  be  ob¬ 
tained  prior  to  the  reconstruction  process.  The  2D  breast 
boundary  curves  are  extracted  from  the  individual  PV  images 
with  a  segmentation  method  developed  in  our  laboratory. 
The  3D  breast  surface  is  generated  by  back  projecting  the  2D 
breast  boundary  curves  sequentially  followed  by  a  3D  coni¬ 
cal  trimming  process.  The  computational  effort  involved  is 
trivial  compared  to  the  reconstruction  procedure  because 


only  boundary  pixels  are  needed  to  be  back  projected.  Note 
that  the  3D  conical  trimming  method  cannot  recover  parts  of 
the  breast  surface  where  they  overlap  projected  breast  re¬ 
gions  in  all  PV  images.  This  will  occur  most  likely  in  the 
region  close  to  the  gap  between  the  breast  support  plate  and 
the  lower  part,  and  between  the  compression  plate  and  the 
upper  part,  around  the  boundary  of  the  compressed  breast,  as 
indicated  by  dashed  arrows  in  Fig.  3.  The  size  of  the  breast 
surface  area  that  cannot  be  recovered  due  to  missing  infor¬ 
mation  will  depend  on  the  angular  range  of  the  DTM  system, 
and  the  thickness  and  boundary  curvature  of  the  compressed 
breast.  The  larger  the  angular  range  of  the  DTM  system,  the 
smaller  the  missing  gap. 

The  efficiency  ratio  is  defined  by  dividing  the  average 
background  area  outside  the  detected  breast  boundary  over 
all  PV  images  by  the  detector  area.  Therefore  the  saving  on 
the  computational  effort  depends  on  the  breast  size  of  indi¬ 
vidual  patient  case  and  is  greater  for  smaller  breasts.  To  es¬ 
timate  the  average  saving  on  computational  effort  in  a  ran¬ 
dom  sample  of  patient  cases,  we  have  calculated  the 
efficiency  ratio  as  the  background  area  to  detector  area  ratio 
for  a  data  set  of  96  two-view  full  field  digital  mammograms 
(FFDMs)  acquired  with  a  GE  Senographe  2000D  system,  the 
detector  of  which  has  the  same  dimension  as  that  of  our 
prototype  DTM  system.  This  data  set  was  collected  previ¬ 
ously  in  the  Department  of  Radiology  at  the  University  of 
Michigan  for  our  CAD  study.4'1  The  resulting  mean  effi¬ 
ciency  ratio  for  96  CC-view  FFDMs  is  65.2%  with  a  stan¬ 
dard  deviation  of  14.7%;  and  the  mean  value  for  96  MLO- 
view  FFDMs  is  56.0%  with  a  standard  deviation  of  13.5%. 
While  the  average  background  area  in  the  21  PVs  of  a  DTM 
will  be  slightly  different  from  that  of  an  FFDM,  this  should 
provide  a  good  estimate  of  the  efficiency  ratio. 

The  developed  2D  and  3D  breast  boundary  information 
can  be  used  in  different  ways  with  DTM  reconstruction  al¬ 
gorithms,  rather  than  the  one  we  described  in  this  work.  For 
example,  with  SART  or  other  iterative  methods,  the  3D 
breast  surface  can  be  employed  at  different  stages  in  the 
iterative  process,  e.g.,  masking  the  resulting  imaged  volume 
to  allow  setting  all  voxels  outside  the  breast  surface  to  air 
after  each  PV  image  has  been  used,  rather  than  after  all  PV 
images  have  been  used  exactly  once,  as  we  suggested  in  this 
work.  In  terms  of  using  object  shape  information  in  image 
reconstruction,  Tam  has  presented  a  simple  method  to  con¬ 
struct  the  convex  hull  of  an  imaged  object  in  2D  limited- 
angle  computerized  tomography  (CT)  by  intersecting  or  su¬ 
perimposing  the  back-projected  object  shadows.47  He 
utilized  this  object  support  information  in  an  iterative 
limited-angle  CT  reconstruction  in  which  the  missing  PV 
images  are  estimated  and  used  in  a  full-range  FBP  method. 
During  the  iterative  process,  the  reconstructed  object  density 

distribution  is  improved  by  correcting  the  external  volume  to 

48 

zero  and  applying  other  constraints  to  the  density  range. 

Another  potential  application  of  the  3D  breast  surface  in¬ 
formation  is  quantitative  DTM  reconstruction.  In  DTM  re¬ 
construction,  the  air  volume  outside  the  breast,  if  not  ex¬ 
cluded,  will  affect  the  estimation  of  the  average  attenuation 
value  back  projected  to  the  real  breast  volume.  With  the  3D 
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breast  surface  information,  if  only  the  pathlength  within  the 
real  breast  volume  will  be  taken  into  account  in  calculating 
the  x-ray  attenuation  properties  of  the  breast  tissues,  it  will 
provide  potential  improvement  in  quantitative  estimation  of 
the  linear  attenuation  coefficients. 

In  addition  to  applications  in  DTM  reconstruction,  the  3D 
breast  surface  and  volume  information  will  be  useful  for 
breast  density  estimation.  Breast  density  has  been  shown  to 
be  an  important  risk  factor  for  breast  cancer.  The  change  in 
breast  density  is  considered  a  useful  surrogate  for  estimating 
the  change  in  breast  cancer  risk  due  to  treatment  or  interven¬ 
tion.  Currently  the  change  in  breast  density  is  monitored 
mainly  by  estimating  the  change  in  the  dense  area  on  con¬ 
ventional  projection  mammograms.  Because  of  the  projec¬ 
tion  of  the  3D  volume  to  a  2D  plane,  the  percentage  of  dense 
area  on  a  mammogram  does  not  accurately  reflect  the  per¬ 
centage  of  dense  tissue  in  the  breast  volume  although  a 
strong  correlation  has  been  demonstrated.49  It  can  be  ex¬ 
pected  that  3D  volumetric  information  available  from  DTM 
will  provide  a  better  estimation  of  the  amount  of  dense  tissue 
in  the  breast  volume  and  its  changes  over  time  than  projec¬ 
tion  mammograms. 

V.  CONCLUSION 

In  this  study,  we  have  applied  the  2D  breast  boundary 
information  and  the  generated  3D  breast  surface  to  SART 
reconstruction  in  breast  tomosynthesis  mammography.  The 
2D  breast  boundary  is  detected  on  the  projection  view  im¬ 
ages  using  a  boundary  tracking  algorithm  developed  in  our 
laboratory  and  the  3D  breast  surface  is  generated  with  a  3D 
conical  trimming  method.  The  2D  breast  boundary  curves  on 
PV  images  are  used  to  restrict  the  SART  reconstruction  to  be 
performed  only  inside  the  breast  area  while  the  3D  breast 
surface  is  used  to  exclude  reconstruction  artifacts  outside  the 
breast  volume.  Experimental  results  with  patient  PV  images 
demonstrated  that  the  proposed  method  can  substantially  im¬ 
prove  computational  efficiency  by  eliminating  unnecessary 
reconstruction  in  regions  outside  the  breast.  Both  breast 
boundary  and  detector  boundary  artifacts  can  be  effectively 
removed  by  the  proposed  method. 
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ABSTRACT 

Digital  Tomosynthesis  Mammography  (DTM)  is  an  emerging  technique  that  has  the  potential  to 
improve  breast  cancer  detection.  DTM  acquires  low-dose  mammograms  at  a  number  of  projection  angles 
over  a  limited  angular  range  and  reconstructs  the  3D  breast  volume.  Due  to  the  limited  number  of  projections 
within  a  limited  angular  range  and  the  finite  size  of  the  detector,  DTM  reconstruction  contains  boundary  and 
truncation  artifacts  that  degrade  the  image  quality  of  the  tomosynthesized  slices,  especially  that  of  the 
boundary  and  truncated  regions.  In  this  work,  we  developed  artifact  reduction  methods  that  make  use  of  both 
2D  and  3D  breast  boundary  information  and  local  intensity-equalization  and  tissue-compensation  techniques. 
A  breast  phantom  containing  test  objects  and  a  selected  DTM  patient  case  were  used  to  evaluate  the  effects 
of  artifact  reduction.  The  contrast-to-noise  ratio  (CNR),  the  normalized  profiles  of  test  objects,  and  a  non¬ 
uniformity  error  index  were  used  as  performance  measures.  A  GE  prototype  DTM  system  was  used  to 
acquire  21  PVs  in  3°  increments  over  a  ±30°  angular  range.  The  Simultaneous  Algebraic  Reconstruction 
Technique  (SART)  was  used  for  DTM  reconstruction.  Our  results  demonstrated  that  the  proposed  methods 
can  improve  the  image  quality  both  qualitatively  and  quantitatively,  resulting  in  increased  CNR  value, 
background  uniformity  and  an  overall  reconstruction  quality  comparable  to  that  without  truncation.  For  the 
selected  DTM  patient  case,  the  obscured  breast  structural  information  near  the  truncated  regions  was 
essentially  recovered.  In  addition,  restricting  SART  reconstruction  to  be  performed  within  the  estimated  3D 
breast  volume  increased  the  computation  efficiency. 

Keywords:  Digital  Tomosynthesis  Mammography  (DTM),  Simultaneous  Algebraic  Reconstruction 
Technique  (SART),  truncation  artifacts,  boundary  artifacts. 

1.  INTRODUCTION 

Digital  tomosynthesis  mammography  (DTM)  is  an  emerging  technique  that  can  provide  volumetric 
information  for  breast  imaging  by  reconstructing  the  breast  volume  from  multiple  low-dose  projection  views 
acquired  at  a  number  of  angles  in  a  limited  angular  range1,2.  The  total  radiation  dose  is  set  to  be  comparable 
to  that  used  in  regular  screening  mammography.  DTM  reconstruction  is  therefore  a  limited-angle  cone-beam 
tomographic  problem.  We  previously  demonstrated  that  the  Simultaneous  Algebraic  Reconstruction 
Technique  (SART)  method  can  achieve  high  image  quality3.  However,  the  reconstructed  DTM  slices  contain 
strong  boundary  and  truncation  artifacts  due  to  the  limited  angular  range  and  the  truncated  projection-view 
(PV)  images.  Moreover,  iterative  reconstruction  algorithm  such  as  SART  requires  extensive  computation, 
therefore  it  is  of  practical  importance  to  improve  computational  efficiency. 
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In  this  work,  our  purpose  is  to  develop  methods  to  reduce  both  truncation  and  boundary  artifacts, 
improve  reconstruction  image  quality  in  the  truncated  region,  and  enhance  computational  efficiency  of  the 
SART  method. 


2.  METHODS 


2.1  DTM  SYSTEM 

A  GE  prototype  DTM  system  was  used  to  acquire  21  PVs  in  3°  increments  over  a  ±30°  angular  range. 
The  imaging  geometry  of  this  DTM  system  is  schematically  shown  in  Figure  1.  Institutional  Review  Board 
(IRB)  approval  and  informed  consent  were  obtained  for  patient  cases.  The  distance  from  the  x-ray  focal  spot 
to  the  center  of  rotation  is  64  cm  and  the  plane  along  which  the  x-ray  source  rotates  is  perpendicular  to  the 
detector  surface  at  the  chest  wall  edge.  The  focal-spot-to-detector  distance  is  66  cm.  The  system  has  a  Csl 
phosphor/a:Si  active  matrix  flat  panel  digital  detector  with  a  pixel  size  of  0.1  mmxQ.l  mm  and  16  bits  per 
pixel.  The  image  acquisition  process  takes  less  than  8  seconds  and  the  digital  detector  (X-Y  plane)  is 
stationary  during  image  acquisition. 

In  our  previous  work3,  we  have  compared  different  methods  for  DTM  reconstruction  based  on  breast 
phantom  study  and  demonstrated  that  the  SART  method5  can  provide  high-quality  tomosynthesized  slices 
while  being  more  efficient  than  the  maximum  likelihood  method2.  In  this  study,  SART  was  employed  for 
DTM  reconstruction.  The  imaged  volume  was  subdivided  into  a  set  of  voxels,  of  which  the  X  and  Y 
dimensions  were  set  to  be  the  same  as  the  detector  pixel  size,  i.e.,  0.1  mm,  and  the  Z  dimension  (slice 
thickness)  was  chosen  to  be  1.0  mm.  In  SART,  an  update  of  the  attenuation  properties  of  the  imaged  volume 
is  performed  based  on  each  PV  and  the  iteration  process  starts  with  an  initial  distribution.  The  value  of  each 
voxel  is  updated  after  all  rays  in  one  projection  view  have  been  processed.  We  used  SART  with  one  iteration, 
i.e.  using  all  PV  images  exactly  once,  in  this  study.  The  relaxation  parameter  was  set  to  be  0.5  for  all  PVs. 

2.2  BOUNDARY  ARTIFACT  REDUCTION 

The  2D  breast  boundary  on  each  PV  image  is  segmented  by  a  breast  boundary  detection  program 
developed  in  our  laboratory6.  There  are  two  main  steps  in  the  algorithm:  detection  of  an  initial  breast 
boundary  by  applying  Otsu’s  thresholding7  to  the  histogram  of  the  PV  image,  followed  by  a  boundary 
tracking  procedure  based  on  gradient  information  obtained  from  both  horizontal  and  vertical  Sobel  filtering. 
With  the  detected  2D  breast  boundary  curves  on  all  PV  images,  a  3D  breast  surface  can  be  generated  to 
enclose  the  compressed  breast  while  excluding  the  exterior  air  space.  We  implemented  a  “3D  conical 
trimming”  process  to  recover  the  unique  convex  hull  inside  the  3D  imaged  volume  of  which  the  projections 
at  different  angles  precisely  correspond  to  the  breast  shadow  on  each  PV  image.  Specifically,  for  a  given  PV 
image,  the  breast  boundary  curve  is  back-projected  into  the  imaged  volume  and  all  “air- volume”  outside  the 
generated  conical  surface  formed  by  the  x-ray  source  and  the  breast  boundary  is  trimmed  off.  This  process  is 
performed  for  all  PV  images.  The  intersection  of  the  conical  surfaces  from  all  PV  breast  boundary  curves 
forms  the  desired  convex  hull.  The  convex  hull  combined  with  the  top  and  bottom  surfaces  of  the  breast 
delineated  by  the  compression  paddle  and  the  breast  support  plate  define  an  enclosed  breast  volume. 

To  apply  the  breast  boundary  information  to  DTM  reconstruction  using  the  SART  method,  we 
restrict  the  re -projection  and  the  back-projection  of  the  difference  between  the  calculated  data  and  the 
measured  data  to  be  performed  along  those  “valid”  rays  within  the  2D  breast  area  only.  After  one  complete 
SART  iteration  has  been  performed,  we  use  the  3D  breast  surface  to  mask  the  resulting  volume  and  set  the 
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values  of  all  voxels  outside  the  breast  surface  to  that  of  air.  This  will  remove  all  boundary  artifacts,  including 
those  from  the  breast  boundary  and  the  detector  boundary  outside  the  breast  volume.  The  computational 
effort  is  reduced  by  eliminating  the  processing  of  all  rays  outside  the  breast  region  on  each  individual  PV. 
Detailed  description  of  this  method  can  be  found  in  our  recent  publication8. 

2.3  TRUNCATION  ARTIFACT  REDUCTION 

When  a  portion  of  the  breast  is  not  recorded  in  some  or  all  of  the  PV  images  because  of  a  finite-size 
detector,  the  truncated  projections  will  cause  artifacts  in  DTM  reconstruction.  The  limited  field  of  view 
(FOV)  results  in  unexposed  regions  in  the  imaged  volume  and  truncation  of  the  PVs,  particularly  at  large 
projection  angles.  In  iterative  reconstruction  such  as  SART,  the  imaged  volume  is  updated  by  processing 
each  individual  PV  image,  i.e.  only  the  part  of  the  imaged  volume  within  the  cone-beam  ray  path  of  the 
current  PV  will  be  updated.  This  will  result  in  discontinuity  in  the  voxel  values  across  the  cone-beam  path 
boundary  which  will  be  enhanced  by  further  PV  processing  in  the  same  and  subsequent  iterations.  These 
artifacts,  referred  to  as  “detector  boundary”  truncation  artifacts,  appear  as  bright  staircase-like  lines  or  bands 
at  the  two  sides  perpendicular  to  the  x-ray  source  motion  on  all  tomosynthesized  slices. 

A  second  source  of  truncation  artifacts  is  caused  by  the  missed  portion  of  breast  tissue  that  is  outside 
the  finite  imaged  volume  modeled  in  the  reconstruction  algorithm.  The  missed  tissue  volume  will  cause 
estimation  error  in  the  x-ray  attenuation.  Since  the  attenuation  will  be  considered  to  have  occurred  only 
within  the  shortened  pathlength,  the  voxel  values  will  be  overestimated.  The  overestimated  attenuation 
results  in  bright  voxels  and  will  be  referred  to  as  “glaring”  artifact. 

Both  of  these  truncation  artifacts  are  apparent  mainly  at  the  image  boundary  of  DTM  slices.  They 
will  obscure  the  breast  tissue  details  near  the  boundary  of  DTMs,  potentially  affecting  the  accuracy  of  lesion 
detection.  These  artifacts  cannot  be  eliminated  by  using  the  breast  boundary  information  alone.  We  have 
developed  correction  techniques  to  reduce  their  visibility. 

For  truncation  artifact  reduction,  a  local  intensity-equalization  method  was  developed  to  suppress  the 
discontinuity  of  the  reconstructed  voxel  intensity.  Specifically,  for  each  PV  image,  after  the  updating  using 
the  current  PV  is  completed,  we  replace  all  voxel  values  within  the  oblique  wedge-shape  volume  with  some 
average  values  obtained  from  their  neighborhood  region.  This  procedure  is  implemented  in  a  column-wise 
manner,  i.e.,  along  each  column  of  voxels  parallel  to  the  Y-axis  on  every  tomosynthesized  slice.  To  reduce 
the  glaring  artifact,  we  designed  a  practical  tissue-compensation  method.  For  any  x-ray  incident  on  the  preset 
rectangular  imaged  volume  from  the  boundary  side,  it  is  assumed  that  this  ray  has  encountered  very  large 
attenuation  if  the  corresponding  detected  x-ray  intensity  is  below  a  predetermined  threshold.  The  computed 
total  attenuation  of  this  ray  will  be  compensated  for  by  assuming  that  the  missing  volume  is  filled  with 
“average”  breast  tissue.  In  addition,  the  modified  difference  between  the  detected  data  and  the  new  computed 
data  is  normalized  by  the  compensated  total  pathlength  but  back-projected  only  to  those  voxels  within  the 
preset  imaged  volume. 

2.4  BREAST  PHANTOM  AND  FIGURES  OF  MERIT 

To  quantitatively  evaluate  the  improvement  in  image  quality  by  artifact  correction,  a  custom-built 
breast  phantom  was  used  which  contained  test  objects  on  a  thin  feature  layer  sandwiched  between 
homogeneous  Lucite  plates.  The  test  objects  were  arranged  in  groups  in  both  the  boundary  and  middle 
regions,  as  schematically  shown  in  Figure  2.  The  test  objects  included  layers  of  thin  circular  aluminum  foils 
to  simulate  masses,  high-contrast  metal  wires  to  simulate  metal  markers,  and  low-contrast  plastic  wires  to 
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simulate  speculations.  The  two  groups  of  masses,  denoted  by  Ml,  M2  and  M3  and  by  M4,  M5  and  M6 
respectively,  differ  in  the  number  of  aluminum  foils  to  simulate  different  contrasts.  The  lower-contrast 
masses  are  drawn  with  a  slightly  smaller  size  than  those  for  the  higher-contrast  ones  although  they  have  the 
same  diameter  in  reality. 

The  contrast-to-noise  ratio  (CNR),  the  normalized  line  profiles  of  feature  objects,  and  a  non¬ 
uniformity  error  index  are  used  as  performance  measures  to  compare  the  image  quality  of  the  test  objects 
near  the  detector  boundary  with  and  without  applying  the  artifact  reduction  methods,  and  to  the  image  quality 
of  the  test  objects  in  the  middle  region  where  no  truncation  occurred.  The  CNR  of  the  simulated  mass  is 
defined  as  the  difference  in  the  average  pixel  intensities  between  the  mass  and  its  local  image  background, 
divided  by  the  standard  deviation  of  pixel  intensity  in  the  same  background.  The  non-uniformity  error  index 
(NUEI)  is  defined  as  the  ratio  of  the  difference  between  the  average  pixel  intensities  of  two  ROIs  located  in 
the  middle  region  and  the  lower  boundary  region,  respectively,  to  the  average  pixel  intensity  of  the  middle 
region.  Both  ROIs  in  the  boundary  and  the  middle  regions  were  chosen  to  be  homogeneous  Lucite 
background  containing  no  objects. 


3.  RESULTS 

For  the  breast  phantom  study,  the  number  of  PV  images  in  which  the  group  of  test  objects  in  the 
lower  boundary  region  (i.e.,  Ml,  Wl,  M4  and  W4)  was  truncated  was  8  to  9.  Figure  3  shows  the 
reconstruction  results  of  the  lower  part  of  the  feature  layer  with  and  without  the  truncation  artifact  reduction. 
Results  of  the  upper  and  middle  groups  of  features  were  not  shown  because  they  were  not  affected  by  the 
proposed  methods.  The  detector  boundary  artifacts  appearing  as  bright  horizontal  lines  have  been  totally 
removed.  The  visual  image  quality  of  all  features  was  markedly  improved,  especially  for  the  low-contrast 
plastic  wire  that  was  almost  totally  obscured  by  the  truncation  artifacts  without  correction.  The  line  profiles 
of  all  simulated  masses  are  shown  in  Figure  4,  in  which  the  high-  and  low-contrast  mass  groups  are  shown 
separately.  In  each  group,  the  line  profiles  of  the  object  located  in  the  lower  region  (e.g.,  Ml),  without  and 
with  artifact  reduction,  were  plotted  together  with  the  line  profiles  of  the  counterpart  object  in  the  middle 
region  (e.g.,  M2)  and  the  top  region  (e.g.,  M3).  For  masses,  all  line  profiles  were  mean-removed.  It  can  be 
seen  that  the  line  profiles  of  all  masses  located  in  the  lower  region  contained  strong  truncation  artifacts  (i.e. 
large  discontinuities)  without  artifact  reduction.  With  artifact  reduction,  the  resulting  line  profiles  were  much 
smoother,  similar  to  that  of  the  counterpart  object  in  the  middle  region  where  no  truncation  occurred.  The 
CNR  values  of  the  simulated  masses  and  the  NUEI  values  are  listed  in  Table  1  where  the  numbers  in  the 
parentheses  are  the  results  without  artifact  reduction.  The  NUEI  values  indicate  that  the  artifact  reduction 
methods  improved  the  uniformity  in  the  reconstructed  homogeneous  Fucite  background.  The  improvement 
in  the  CNR  values  of  masses  Ml  and  M4  is  substantial. 

For  the  selected  patient  DTM  case,  a  mediolateral  oblique  (MFO)  view  with  truncated  PV  images 
was  used.  The  #1,  #11  and  #21  raw  PV  images,  corresponding  to  the  projection  angles  -30°,  0°,  +30°,  are 
shown  in  Figures  5.  The  x-ray  source  moved  in  the  vertical  direction  relative  to  the  images  shown.  Figure  6 
shows  the  generated  3D  breast  surface  by  using  the  proposed  3D  conical  trimming  method.  Figure  7  shows 
examples  of  DTM  slices  reconstructed  with  the  SART  method  for  which  a  constant  distribution  of  voxel 
values  was  used  as  initialization.  The  DTM  images  reconstructed  without  using  any  artifact  reduction  (top 
row)  and  with  both  boundary  and  truncation  artifact  reduction  (bottom  row)  were  compared.  The  three 
chosen  slices  are  3.1  cm,  2.1  cm,  and  0.8  cm  above  the  breast  support  plate  and  the  thickness  of  the  imaged 
breast  was  measured  by  the  DTM  system  as  4.9  cm.  Without  using  any  artifact  reduction,  the  detector 
boundary  artifacts  were  observed  as  multiple  horizontal  lines  especially  at  the  top  of  all  three 
tomosynthesized  slices.  The  breast  boundary  artifacts,  represented  as  repeated  breast-shape  ghosts  around  the 
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boundary,  were  clearly  observed  around  the  true  breast  boundary  in-focus  on  each  slice.  It  is  also  evident  that 
both  the  detector  boundary  artifact  and  glaring  artifact  strongly  degraded  the  visibility  of  the  breast  structures 
in  the  boundary  areas  that  were  close  to  the  pectoral  muscle  where  the  PV  images  were  severely  truncated. 
The  correction  methods  removed  all  breast  boundary  artifacts,  effectively  suppressed  the  truncation  artifacts 
and  recovered  the  local  breast  tissue  details.  The  glaring  artifacts  caused  by  the  truncated  pathlengths  have 
been  corrected  and  the  reconstructed  voxel  intensities  are  similar  to  those  of  normal  breast  tissue.  In  this 
case,  the  computation  time  was  reduced  by  63.2%  by  using  the  2D  breast  boundary  masks,  compared  to  full- 
field  processing. 


4.  DISCUSSION  AND  CONCLUSION 

DTM  can  provide  quasi-3D  information  in  breast  imaging  in  contrast  to  conventional  mammography. 
However,  due  to  the  limited  number  of  projections  within  a  limited  angular  range  and  the  finite  size  of  the 
flat  panel  detector,  DTM  reconstruction  contains  boundary  and  truncation  artifacts  that  degrade  the  image 
quality  of  the  tomosynthesized  slices,  especially  in  the  boundary  and  truncated  regions.  In  this  study,  we 
developed  artifact  reduction  methods  and  evaluated  the  improvement  on  breast  phantom  images  and  DTM 
patient  case.  Our  results  indicated  that,  based  on  the  SART  method  for  the  limited-angle  cone-beam  3D 
reconstruction,  the  proposed  methods  can  substantially  reduce  the  artifacts,  improve  the  image  quality  both 
qualitatively  and  quantitatively,  and  provide  reconstruction  quality  comparable  to  that  without  truncation.  In 
addition,  the  computation  time  of  SART  is  reduced  substantially  by  using  the  2D  breast  boundary 
information.  Further  work  is  underway  to  investigate  the  generalizability  of  the  correction  methods  to  other 
reconstruction  techniques  and  the  correction  of  artifacts  that  are  caused  by  high-attenuation  objects. 
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Figure  3.  SART  reconstructed  DTM  slice  of  the  phantom  feature  layer,  without  (left)  and  with  (right)  the 
artifact  reduction.  Only  the  lower  region  was  shown  and  the  corresponding  test  objects  are  indicated. 


(a)  Ml,  M2  and  M3 


Pixel  Index  (0.1  mm/pixel) 

(b)  M4,  M5  and  M6 


Figure  4.  Comparison  of  the  normalized  line  profiles  of  all  simulated  masses:  the  feature  located  in  the 
lower  region,  without  (solid  gray  line)  and  with  artifact  reduction  (solid  black  line),  the 
corresponding  feature  located  in  the  middle  (black  dashed  line)  and  the  top  (black  dotted  line) 
regions.  The  pixel  intensity  was  normalized  by  removing  the  mean  of  each  line  profile. 
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Figure  5.  Projection  view  images  of  a  DTM  patient  case  in  MLO  view:  PV  at  -30°  (left),  PV  at  0° 

(middle)  and  PV  at  +30°  (right).  The  x-ray  source  moved  in  the  vertical  direction  relative  to  the 
images  shown. 


Figure  6.  The  3D  breast  surface  generated  from  2D  breast  boundary  curves  by  using  the  3D 
conical  trimming  method. 
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ABSTRACT 

Digital  tomosynthesis  mammography  (DTM)  can  provide  quasi-3D  structural  information  of  the  breast  by 
reconstructing  the  breast  volume  from  projection  views  (PV)  acquired  in  a  limited  angular  range.  Our  purpose  is  to 
design  an  effective  classifier  to  distinguish  breast  masses  from  normal  tissues  in  DTMs.  A  data  set  of  100  DTM  cases 
collected  with  a  GE  first  generation  prototype  DTM  system  at  the  Massachusetts  General  Hospital  was  used.  We 
reconstructed  the  DTMs  using  a  simultaneous  algebraic  reconstruction  technique  (SART).  Mass  candidates  were 
identified  by  3D  gradient  field  analysis.  Three  approaches  to  distinguish  breast  masses  from  normal  tissues  were 
evaluated.  In  the  3D  approach,  we  extracted  morphological  and  run-length  statistics  texture  features  from  DTM  slices  as 
input  to  a  linear  discriminant  analysis  (LDA)  classifier.  In  the  2D  approach,  the  raw  input  PVs  were  first  preprocessed 
with  a  Laplacian  pyramid  multi-resolution  enhancement  scheme.  A  mass  candidate  was  then  forward-projected  to  the 
preprocessed  PVs  in  order  to  determine  the  corresponding  regions  of  interest  (ROIs).  Spatial  gray-level  dependence 
(SGLD)  texture  features  were  extracted  from  each  ROI  and  averaged  over  11  PVs.  An  LDA  classifier  was  designed  to 
distinguish  the  masses  from  normal  tissues.  In  the  combined  approach,  the  LDA  scores  from  the  3D  and  2D  approaches 
were  averaged  to  generate  a  mass  likelihood  score  for  each  candidate.  The  A:  values  were  0.87±0.02,  0.86±0.02,  and 
0.91±0.02  for  the  3D,  2D,  and  combined  approaches,  respectively.  The  difference  between  the  Az  values  of  the  3D  and 
2D  approaches  did  not  achieve  statistical  significance.  The  performance  of  the  combined  approach  was  significantly 
(p<0.05)  better  than  either  the  3D  or  2D  approach  alone.  The  combined  classifier  will  be  useful  for  false-positive 
reduction  in  computerized  mass  detection  in  DTM. 

Keywords:  Digital  Tomosynthesis  Mammography  (DTM),  Simultaneous  Algebraic  Reconstruction  Technique 

(SART),  Breast  mass,  Receiver  operating  characteristic  (ROC)  analysis 


1.  INTRODUCTION 

Breast  cancer  is  one  of  the  leading  causes  of  cancer  mortality  among  women1.  Although  mammography  is  a  powerful 
screening  tool  for  detecting  breast  cancer,  the  sensitivity  of  cancer  detection  is  often  limited  by  the  presence  of  the 
overlapping  dense  fibroglandular  tissue  in  the  breast.  Digital  tomosynthesis  mammography  (DTM)  is  one  of  the  new 
breast  imaging  modalities  that  have  potential  to  improve  the  detection  and  diagnosis  of  breast  cancer.  In  DTM,  a  serial 
of  low-dose  projection  view  images  (PVs)  are  acquired  when  the  x-ray  source  is  rotated  about  the  fulcrum  over  a  limited 
angular  range.  DTM  slices  reconstructed  from  the  serial  of  PVs  can  provide  quasi-3D  structural  information  of  the 
breast  which  may  improve  the  conspicuity  of  breast  cancer  and  also  reduce  the  effects  of  the  overlapping  dense  tissues. 


In  conventional  mammography,  computer-aided  detection  system  (CADe)  can  increase  the  breast  cancer  detection  rate 
by  radiologists  both  in  the  laboratory  and  in  clinical  practice2"7.  We  are  developing  a  CADe  system  to  assist 
radiologists  in  detecting  breast  masses  in  DTMs.  In  this  preliminary  study,  our  purpose  is  to  design  an  effective 
classifier  to  distinguish  breast  masses  from  normal  tissues  for  false  positive  (FP)  reduction  in  computerized  mass 
detection  in  DTMs. 
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Figure  1.  Histogram  of  the  mass  sizes  measured  on  the  central  slices  of  the  reconstructed  DTMs.  Mass  sizes  are 
measured  as  the  longest  dimension  of  the  mass  by  an  experienced  MQSA  radiologist.  The  size  of  the 
masses  in  this  data  set  ranged  from  5.5  mm  to  43.4  mm  with  a  mean  size  of  17.4±7.3  mm. 
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Figure  2.  An  example  of  a  spiculated  mass  in  DTMs.  (a)  PV  at  0  degree,  (b)  the  DTM  slice  approximately  intersecting  the  center 
of  the  mass,  (c)  ROI  image  on  PV  at  0  degree,  and  (d)  ROI  image  on  the  DTM  slice  in  (b). 
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2.  MATERIALS  AND  METHODS 


2.1  Materials 

DTMs  in  this  study  were  collected  with  Institutional  Review  Board  (IRB)  approval  in  the  Breast  Imaging  Research 
Laboratory  at  the  Massachusetts  General  Hospital  with  a  GE  first  generation  prototype  DTM  system.  The  DTM  system 
has  a  flat  panel  detector  with  a  pixel  size  of  0.1  mm  x  0.1  mm.  It  acquired  11  PVs  in  5-degree  increments  over  a  total 
tomosynthesis  angle  of  50-degree.  The  total  dose  for  the  1 1  PVs  was  designed  to  be  less  than  1.5  times  that  of  a  single 
conventional  mammogram.  Only  mediolateral  oblique  view  (MLO)  DTM  is  collected  for  each  case.  The  data  set  we 
used  in  this  study  contained  100  cases.  Each  case  had  one  single  breast  DTM  which  was  reconstructed  by  a 
simultaneous  algebraic  reconstruction  technique  (SART)  developed  in  previous  study8.  SART  method  iteratively 
updates  the  imaged  volume  using  each  PV  image  sequentially.  The  number  of  reconstructed  slices  in  this  data  set 
ranged  from  45  to  90  depending  on  the  compressed  breast  thickness.  An  experienced  MQSA  radiologist  identified  the 
location  of  the  masses  in  the  reconstructed  DTM  slices.  There  were  69  malignant  and  3 1  benign  masses  in  this  data  set 
including  some  large-area  architectural  distortions.  Figures  1  showed  the  histograms  of  mass  sizes  measured  on  the 
central  slices.  The  mass  size  ranged  from  5.5  to  43.4  mm  with  a  mean  size  of  17.4±7.3  mm  on  the  central  DTM  slices. 
Figure  2  is  an  example  of  a  DTM  with  a  spiculated  mass;  the  PV  at  0  degree  and  the  SART  reconstructed  slice 
approximately  intersecting  the  center  of  the  mass  are  shown. 

2.2  Methods 

2.2.1  CADe  System  Overview 


Figure  3.  Block  diagram  of  the  CADe  system  for  mass  detection  on  DTMs. 


Our  CADe  system  consists  of  five  processing  steps:  1)  segmentation  of  breast  region  on  each  reconstructed  slice,  2)  pre¬ 
screening  of  mass  candidates,  3)  identification  of  suspicious  objects,  4)  extraction  of  morphological  and  texture  features, 
and  5)  classification  between  the  normal  tissue  and  the  mass  regions  by  using  LDA  classifiers.  The  block  diagram  for 
the  CADe  system  is  shown  in  Figure  3. 
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We  have  described  a  prototype  CADe  system  for  detection  of  masses  in  the  reconstructed  DTM  volumes  previously.9 
For  each  reconstructed  slice,  the  CADe  system  first  segments  the  breast  region  from  the  background  by  using  a  breast 
boundary  detection  algorithm.  The  rest  of  the  processes  is  only  applied  to  the  segmented  breast  region.  For  the  pre¬ 
screening  stage,  we  designed  a  3D  filter  using  the  information  of  gradient  field  directions  to  identify  mass  candidates. 
After  this  enhancement  filtering,  the  local  maxima  within  the  breast  region  were  identified  as  the  mass  candidates  on 
each  DTM  volume.  Volumes  of  interest  (VOI)  were  then  identified  in  the  breast  with  the  center  of  a  VOI  placed  at 
each  of  the  mass  candidates.  The  suspicious  structure  in  each  identified  location  was  extracted  by  a  3D  region  growing 
method. 

During  the  feature  analysis  and  classification  step,  we  evaluated  three  approaches  to  distinguish  breast  masses  from 
normal  tissues.  In  the  first  approach,  referred  to  as  3D  approach,  we  extracted  morphological  and  run-length  statistics 
texture  features  from  the  DTM  slices.  A  linear  discriminant  analysis  (LDA)  classifier  with  stepwise  feature  selection 
was  then  trained  to  merge  the  useful  features.  In  the  second  approach,  referred  to  as  2D  approach,  the  raw  input  PVs 
were  first  preprocessed  with  a  Laplacian  pyramid  multiresolution  enhancement  scheme.  A  mass  candidate  was  then 
forward-projected  onto  the  each  of  the  11  PVs  in  order  to  determine  the  corresponding  regions  of  interest  (ROIs)  on  the 
preprocessed  PVs.  Spatial  gray-level  dependence  (SGLD)  texture  features  were  extracted  from  each  ROI  and  averaged 
over  the  11  PVs.  An  LDA  classifier  with  stepwise  feature  selection  was  designed  to  distinguish  the  masses  from 
normal  tissues.  In  the  third  approach,  referred  to  as  the  combined  approach,  the  LDA  scores  of  a  detected  object,  which 
represented  the  likelihood  of  the  object  being  a  mass,  from  the  first  and  second  approaches  were  merged  to  generate  a 
combined  mass  likelihood  score  for  each  mass  candidate.  In  this  study,  a  simple  averaging  method  was  used  to  merge 
the  information. 


2.2.2  Training  and  test  CADe  system 

The  leave -one-case  out  re-sampling  method  was  used  for  training  and  testing  the  FP  classifiers.  During  training, 
feature  selection  with  stepwise  LDA  was  employed  to  obtain  the  best  feature  subset  and  reduce  the  dimensionality  of  the 
feature  space  to  design  an  effective  classifier.  Briefly,  for  a  data  set  of  n  cases,  all  mass  candidates  from  a  test  case  are 
left  out  while  the  other  (n-1)  cases  are  used  for  selection  of  predictor  variables  from  the  feature  pool  and  estimation  of 
the  LDA  classifier  weights  in  each  leave-one-case-out  cycle.  Since  the  appropriate  threshold  values  for  feature  entry, 
feature  elimination,  and  tolerance  of  feature  correlation  were  unknown,  we  used  an  automated  simplex  optimization 
method  to  search  for  the  best  combination  of  thresholds  in  the  parameter  space.  The  trained  LDA  classifier  was  then 
applied  to  the  mass  candidates  in  the  left-out  case  to  obtain  the  mass  likelihood  score  of  each  candidate.  This  process 
was  performed  with  each  of  the  n  cases  left  out  in  turn  so  that  all  objects  would  be  assigned  a  mass  likelihood  score  at 
the  completion  of  the  n  cycles. 


2.2.3  Evaluation  methods 

In  order  to  compare  the  performance  of  the  three  feature  analysis  and  classification  methods,  we  employed  the 
receiver  operating  characteristic  (ROC)  analysis.  The  ROCKIT  curve  fitting  software  and  statistical  significance  test 
for  ROC  analysis  developed  by  Metz  et  al.10  is  used  in  this  study.  The  discriminant  scores  from  the  three  different 
approaches  were  input  as  the  decision  variable  in  the  ROCKIT  program,  which  fits  a  binormal  ROC  curve  based  on 
maximum  likelihood  estimation. 


3.  EXPERIMENTAL  RESULTS 

Figure  4  shows  the  test  ROC  curves  for  classification  of  breast  masses  and  normal  tissues  on  DTMs.  The  areas  under 
the  ROC  curves  (Az)  were  0.87±0.02,  0.86±0.02  and  0.91±0.02  for  the  3D,  2D,  and  combined  approaches,  respectively. 
The  difference  between  the  A-  values  of  the  3D  and  2D  approaches  did  not  achieve  statistical  significance  (p=0.50 ). 
The  performance  of  the  combined  approach  was  significantly  better  than  either  the  3D  or  the  2D  approach  (p=0.04  and 
p=0.03,  respectively).  Table  1  summarized  the  ROC  results. 
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False  Positive  Fraction 

Figure  4.  The  test  ROC  curves  for  classification  of  breast  masses  from  nonnal  tissues  on  digital  tomosynthesis 
mammography  (DTM).  The  A-  values  are  0.87±0.02,  0.86±0.02  and  0.91±0.02  for  the  3D,  2D  and  combined 
approaches,  respectively. 


Table  1.  Estimation  of  the  statistical  significance  in  the  difference  between  the  ROC  performances  of  three  feature 
analysis  and  classification  approaches.  A:  3D  approach,  B:  2D  approach,  C:  combined  approach. 


Az  Value 

P  Value 

A 

0.87 

A  vs  B 

=  0.5 

B 

0.86 

B  vs  C 

<0.05 

C 

0.91 

C  vs  A 

<0.05 

4.  DISCUSSION  AND  CONCLUSIONS 

In  mammography,  3D  anatomical  structures  are  projected  onto  a  2D  image  plane,  and  the  overlapping  ftbroglandular 
breast  tissues  reduce  the  visibility  of  breast  cancer.  DTM  is  one  of  the  promising  methods  that  provide  quasi-3D 
structural  information  and  may  improve  the  sensitivity  and  specificity  of  breast  cancer  detection.  In  this  study,  we 
compared  three  methods  for  classification  of  breast  masses  and  normal  tissues  in  DTMs.  The  forward  projection  for 
extraction  of  2D  image  information  from  the  PVs  is  a  new  method  for  identifying  corresponding  ROIs  for  the  detected 
objects.  The  combined  approach  which  fuses  the  information  on  both  the  2D  projection  views  and  3D  reconstructed 
slices  is  a  promising  approach  to  improving  the  differentiation  of  breast  masses  from  normal  tissues  on  DTMs.  It  is 
expected  that  the  classifier  will  be  useful  for  FP  reduction  in  a  computerized  mass  detection  system. 
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Digital  tomosynthesis  mammography:  Improvement  of  artifact  reduction 
method  for  high- attenuation  objects  on  reconstructed  slices 
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ABSTRACT 

One  major  image  quality  problem  in  digital  tomosynthesis  mammography  (DTM)  is  the  poor  depth-resolution 
caused  by  the  inherent  incomplete  sampling.  This  problem  is  more  pronounced  if  high-attenuation  objects,  such  as 
metallic  markers  and  dense  calcifications,  are  present  in  the  breast.  Strong  ghosting  artifacts  will  be  generated  in  the 
depth  direction  in  the  reconstructed  volume.  Incomplete  sampling  of  DTM  can  also  cause  visible  ghosting  artifacts  in 
the  x-ray  source  motion  direction  on  the  off-focus  planes  of  the  objects.  These  artifacts  may  interfere  with  radiologists’ 
visual  assessment  and  computerized  analysis  of  subtle  mammographic  features.  We  previously  developed  an  artifact 
reduction  method  by  using  3D  geometrical  information  of  the  objects  estimated  from  the  reconstructed  slices.  In  this 
study,  we  examined  the  effect  of  imaging  system  blur  in  DTM  caused  by  the  focal  spot  and  the  detector  modulation 
transformation  function  (MTF).  The  focal  spot  was  simulated  as  a  0.3  mm  square  array  of  x-ray  point  sources.  The 
detector  MTF  was  simulated  using  the  Burgess  model  with  parameters  derived  from  published  data  of  a  GE  FFDM 
detector.  The  spatial-variant  impulse  responses  for  the  DTM  imaging  system,  which  are  required  in  our  artifact 
reduction  method,  were  then  computed  from  the  DTM  imaging  model  and  a  given  reconstruction  technique.  Our  results 
demonstrated  that  inclusion  of  imaging  system  blur  improved  the  performance  of  our  artifact  reduction  method  in  terms 
of  the  visual  quality  of  the  corrected  objects.  The  detector  MTF  had  stronger  effects  than  focal  spot  blur  on  artifact 
reduction  under  the  imaging  geometry  used.  Further  work  is  underway  to  investigate  the  effects  from  other  DTM 
imaging  parameters,  such  as  x-ray  scattering,  different  polyenergetic  x-ray  spectra,  and  different  configurations  of 
angular  range  and  angular  sampling  interval. 

Keywords:  Digital  tomosynthesis  mammography,  tomographic  reconstruction,  artifact  reduction,  constrained  iterative 
method,  imaging  system  blur 


1.  INTRODUCTION 

Digital  tomosynthesis  mammography1'"  (DTM)  capitalizes  on  the  strengths  and  proven  abilities  of  conventional 
mammography,  while  providing  the  much-needed  three-dimensional  (3D)  information.  Display  of  the  tomosynthesized 
breast  volume  with  each  reconstructed  slice  in  sharp  focus  enhances  lesion  visibility  and  facilitates  analysis  of  lesion 
margins  by  reducing  overlapping  structures  in  the  breast  volume.  This  may  allow  the  radiologist  to  not  only  reduce 
missed  cancers,  but  also  make  more  accurate  diagnoses  of  breast  lesions.  One  major  image  quality  problem  in  DTM  is 
the  poor  depth  resolution  caused  by  the  inherent  incomplete  sampling.  This  problem  is  more  pronounced  if  high- 
attenuation  objects,  such  as  metallic  markers  and  dense  calcifications,  are  present  in  the  breast.  Strong  ghosting  artifacts 
will  be  generated  in  the  depth  direction  in  the  reconstructed  volume.  Incomplete  sampling  of  DTM  can  also  cause 
visible  ghosting  artifacts  in  the  x-ray  source  motion  direction  on  the  off-focus  planes  of  the  objects.  These  artifacts  may 
interfere  with  radiologists’  visual  assessment  and  computerized  analysis  of  subtle  mammographic  features. 

A  number  of  methods3'*  have  been  proposed  to  suppress  artifacts  and  enhance  image  quality  for  tomosynthesis. 
We  previously  developed  a  ghosting  artifact  reduction  method  by  using  3D  geometrical  information  of  the  objects 
estimated  from  the  reconstructed  slices9.  The  idea  is  based  on  the  observation  that  high-attenuation  objects  in  breasts 
being  imaged  are  typically  embedded  in  a  local  soft-tissue  background  that  can  be  assumed  to  be  quasi-homogeneous 
within  the  small  local  volume  of  interest  (VOI).  A  set  of  linear  equations  can  be  used  to  describe  the  imaging  process 
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map.  The  prescreening  process  of  an  automated  calcification  detection  system10  for  full-field  digital  mammograms 
(FFDMs)  previously  developed  in  our  laboratory  was  used  for  the  purpose.  In  order  to  find  the  3D  location  and  an 
approximate  volume  of  the  high-attenuation  objects  from  the  2D  segmentation  on  PVs,  a  simple  voting  scheme  was  used. 
In  this  scheme,  each  of  the  voxels  in  the  3D  volume  (0. 1  mm  isotropic  voxel  size)  obtains  one  vote  from  a  PV  if  a  ray 
connecting  the  source  and  a  non-zero  pixel  on  the  2D  binary  map  for  this  PV  passes  through  the  voxel.  Thus  each  of  the 
voxels  in  the  3D  volume  contained  the  number  of  the  votes  out  of  the  21  PV  binary  maps.  A  threshold  of  14  was 
empirically  chosen  for  this  voting  scheme.  The  connected  voxels  with  values  above  14  were  segmented  and  identified  as 
belonging  to  the  same  object.  The  center-of-mass  of  each  segmented  3D  object  was  defined  as  the  object  center.  The 
bounding  box  for  each  3D  object  was  extended  by  3  pixels  (0.3  mm)  on  each  side  on  the  x-y  plane  to  define  the  volume 
of  interest  (VOI)  Q  for  this  object.  If  the  neighboring  VOIs  overlapped,  the  two  objects  would  be  considered  as  one 
and  the  VOI  would  be  extended  to  enclose  both  objects.  Figure  3  shows  the  flowchart  of  the  steps  in  this  stage. 


Fig.  3.  Flowchart  of  the  estimation  of  3D  geometrical  information  of  the  high-attenuation  objects. 

At  the  second  stage,  the  VOI  Q  containing  each  high-attenuation  object  is  modeled  as  a  homogeneous 
background  with  a  high-attenuation  object.  The  eLAC  for  a  voxel  within  Q  is  denoted  as  x  ,p  e  Q  •  The  eLAC  is  the 

sum  of  two  parts,  namely,  the  slow-varying  background  Vp  and  the  contrast  dp-  We  assume  the  VOI  Q  is  isolated  in 

the  sense  that  the  artifacts  outside  Q  caused  by  the  high  eLACs  within  Q  are  small  enough  to  be  ignored.  The 
projection-reconstruction  (PR)  system  within  Q  is  then  characterized  by  a  set  of  linear  equations x  -  x  =  C,,d ,  where  x 

is  the  vector  of  the  voxel  values  %  at  p  within  Q  on  the  reconstruction  slices,  xn  is  the  average  of  the  neighboring 

voxel  values  of  the  high-attenuation  object  to  approximate  the  slow-varying  part  of  xp ,  C_,,  is  the  shift-variant  3D 

impulse  response  matrix  for  the  PR  system  within  Q ,  and  d  is  the  vector  of  dp ’s.  Since  C„  is  generally  sparse  and 

singular  for  DTM,  direct  matrix  inversion  is  time-consuming  and  sensitive  to  the  noise.  Therefore,  a  constrained 
iterative  method11  was  used  to  estimate  the  contrast  eLAC  value  dp  s.  The  3D  segmentation  of  the  high-attenuation 

object  obtained  at  the  previous  stage  was  used  as  the  initialization  in  this  iterative  method.  The  positivity  of  the  eLACs 
was  used  as  the  constraint.  Finally,  the  interplane  artifacts  are  calculated  as  a  linear  combination  of  the  impulse 
responses  corresponding  to  each  of  the  voxels  of  the  high-attenuation  objects  and  subtracted  from  the  original 
reconstructed  slices.  For  the  voxels  associated  with  the  high-attenuation  objects,  the  voxel  values  were  replaced  with  the 
sum  of  the  estimated  background  part  and  the  contrast  part:  xn  +  dp  ■ 
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3.  MODELING  OF  DTM  IMAGING  SYSTEM  BLUR 


There  are  a  number  of  sources  of  blur  in  X-ray  imaging  systems.  In  this  study,  we  examined  the  effects  of 
imaging  blur  in  DTM  system  caused  by  the  focal  spot  and  the  detector  modulation  transformation  function  (MTF).  The 
real  x-ray  source  is  not  a  perfect  point  source.  We  divided  it  into  elements  that  were  small  enough  such  that  we  can 
approximate  its  response  by  a  point.  We  modeled  the  focal  spot  as  a  0.3  mm  array  (3  by  3)  of  point  x-ray  sources  to 
simulate  the  focal  spot  blur.  The  photon  fluence  was  distributed  to  each  simulated  point  x-ray  source  uniformly.  The 
detector  MTF  was  simulated  using  the  Burgess  model12  with  parameters  derived  from  published  data  of  a  GE  FFDM 
detector13.  The  Burgess  model  is  given  by  MTF(u)  =  0.5eifc[aln(u/ u0)\,  where  u  is  the  spatial  frequency  and  erfc[-\  is 

the  complementary  error  function12.  The  two  adjustable  parameters  a  and  u0  were  estimated  as  0.71  and  2.92, 

respectively,  by  curve  fitting.  The  MTF  data  of  the  GE  FFDM  detector13  and  the  MTF  curve  fitted  with  Burgess  Model12 
are  shown  in  Fig.  4.  The  fitted  MTF  is  then  applied  to  the  calculation  of  the  shift-variant  3D  impulse  response  matrix 
Ca  for  the  PR  system  within  Q . 


0.0  2.5  5.0  7.5  10.0  12.5  15.0 


Spatial  frequency  (mm'1) 

Fig.  4.  MTF  of  GE  Senographe  DS  system  (dashed  line)  and  the  MTF  fitted  with  Burgess  Model  (solid  line). 


4.  RESULTS 

A  patient  DTM  with  a  large  dense  calcification  was  selected  from  our  collected  DTM  database  to  evaluate  the 
performance  of  artifact  reduction.  A  simulated  high-attenuation  object  with  preset  eLACs  that  approximately  match  the 
attenuation  of  calcifications  were  digitally  generated  and  projected  onto  the  PVs  near  the  real  dense  calcification  of  the 
patient  DTM.  DTM  slices  were  then  reconstructed  from  the  PVs  and  processed  with  our  artifact  reduction  method  as 
described  in  Section  3.  The  artifact  reduction  performance  for  the  simulated  high-attenuation  object  was  compared  to 
that  of  the  real  calcification  in  the  artifact-reduced  DTMs.  This  artifact  correction  process  was  performed  both  with  and 
without  taking  into  consideration  the  system  blur. 

The  artifact  reduction  results  in  Fig.  5  shows  that  the  inclusion  of  imaging  system  blur  resulted  in  improved 
artifact-reduced  images  for  both  the  simulated  and  real  objects.  The  corrected  objects  are  not  as  artificially  sharp  and 
bright  as  the  images  obtained  from  the  correction  without  considering  the  system  blur.  Some  internal  features  of  the  real 
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ABSTRACT 

We  are  developing  computer-aided  diagnosis  (CADx)  methods  for  classification  of  masses  on  digital  breast 
tomosynthesis  mammograms  (DBTs).  A  DBT  data  set  containing  107  masses  (56  malignant  and  51  benign)  collected  at 
the  Massachusetts  General  Hospital  was  used.  The  DBTs  were  obtained  with  a  GE  prototype  system  which  acquired  11 
projection  views  (PVs)  over  a  50-degree  arc.  We  reconstructed  the  DBTs  at  1-mm  slice  interval  using  a  simultaneous 
algebraic  reconstruction  technique.  The  regions  of  interest  (ROIs)  containing  the  masses  in  the  DBT  volume  and  the 
corresponding  ROIs  on  the  PVs  were  identified.  The  mass  on  each  slice  or  each  PV  was  segmented  by  an  active  contour 
model.  Spiculation  measures,  texture  features,  and  morphological  features  were  extracted  from  the  segmented  mass. 
Four  feature  spaces  were  formed:  (1)  features  from  the  central  DBT  slice,  (2)  average  features  from  5  DBT  slices 
centered  at  the  central  slice,  (3)  features  from  the  central  PV,  and  (4)  average  features  from  all  11  PVs.  In  each  feature 
space,  a  linear  discriminant  analysis  classifier  with  stepwise  feature  selection  was  trained  and  tested  using  a  two  loop 
leave-one-case-out  procedure.  The  test  Az  of  0.91±0.03  from  the  5-DBT-slice  feature  space  was  significantly  (p=0.003) 
higher  than  that  of  0.84±0.04  from  the  1-DBT-slice  feature  space.  The  test  Az  of  0.83±0.04  from  the  11 -PV  feature 
space  was  not  significantly  different  (p=0.18)  from  that  of  0.79±0.04  from  the  1-PV  feature  space.  The  classification 
accuracy  in  the  5-DBT-slice  feature  space  was  significantly  better  (p=0.006)  than  that  in  the  11-PV  feature  space.  The 
results  demonstrate  that  the  features  of  breast  lesions  extracted  from  the  DBT  slices  may  provide  higher  classification 
accuracy  than  those  from  the  PV  images. 

Keywords:  Digital  breast  tomosynthesis,  computer-aided  classification,  masses,  SART 


1.  INTRODUCTION 

In  mammography,  the  presence  of  overlapping  dense  fibroglandular  tissue  reduces  the  conspicuity  of  the  abnormalities, 
making  it  difficult  to  detect  or  to  characterize  a  lesion.  This  problem  may  be  alleviated  by  digital  breast  tomosynthesis 
mammography  (DBT).  In  DBT,  a  series  of  projection  view  (PV)  images  is  acquired  as  the  x-ray  source  is  rotated  over  a 
limited  range  of  angles.  The  total  dose  required  for  DBT  is  kept  at  nearly  the  same  or  only  slightly  higher  than  that  of  a 
regular  mammogram.  Tomographic  slices  of  the  imaged  volume  can  be  generated  with  reconstruction  techniques  from 
the  series  of  PV  images.  The  DBT  slices  provide  quasi-3D  structural  information  and  may  reduce  the  camouflaging 
effects  of  fibroglandular  tissues.  DBT  is  one  of  the  promising  methods  that  may  improve  the  detection  and 
characterization  for  breast  lesions. 1 

Computer-aided  detection  (CADe)  and  computer-aided  diagnosis  (CADx)  have  been  shown  to  improve  breast  cancer 
detection  and  characterization  in  mammography.  It  is  expected  that  computer-aided  image  analysis  will  play  a  role  in 
DBT  because  of  the  large  number  of  images  that  need  to  be  interpreted.  Efforts  are  underway  to  develop  computerized 
detection  systems  for  DBT.2'8  We  have  previously  investigated  the  feasibility  of  developing  a  CADx  system  for  breast 
masses  using  the  reconstructed  DBT  slices  as  input.9  The  purpose  of  the  current  study  is  to  evaluate  an  alternative 
approach  of  using  the  2D  PV  images  as  input,  and  compare  its  performance  with  that  of  the  approach  of  using  the  DBT 
slices  as  input. 
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2.  MATERIALS  AND  METHODS 


2.1  Data  Set 

A  data  set  of  107  masses  (56  malignant  and  51  benign)  was  used  in  this  study.  The  DBTs  were  acquired  with  a  GE 
prototype  system  at  the  Massachusetts  General  Hospital  with  IRB  approval.  The  DBT  system  acquired  11  PVs  in  5- 
degree  increments  over  an  arc  of  50  degrees.  The  breasts  were  imaged  in  the  mediolateral  oblique  (MLO)  view.  We 
reconstructed  the  DBTs  at  1-mm  slice  interval  using  a  simultaneous  algebraic  reconstruction  technique  (SART).  The 
volume  of  interest  (VOI)  containing  the  mass  was  identified  by  a  Mammography  Quality  Standards  Act  (MQSA)- 
approved  radiologist  on  each  DBT  volume.  The  VOI  was  marked  as  a  rectangular  region  of  interest  (ROI)  on  the 
approximately  central  slice  where  the  mass  was  best  visualized,  and  bounded  from  the  top  and  the  bottom  by  the  first 
and  the  last  slice  where  the  mass  became  almost  invisible.  The  mass  ROI  was  extracted  from  each  DBT  slice  with  a 
margin  of  about  4  mm  around  the  bounding  box.  For  a  given  mass,  the  corresponding  ROIs  on  the  PVs  were  located  by 
forward  projection  of  the  central  slice  of  the  mass  to  the  PVs  using  the  known  geometry  of  the  DBT  system.  The  longest 
diameter  of  the  mass  on  the  central  DBT  slice  was  estimated  by  the  radiologist  using  an  electronic  caliper. 

2.2  Feature  Extraction 

The  automated  mass  classification  scheme  used  for  analysis  of  both  DBT  slices  and  PV  images  in  this  study  is  shown  in 
Fig.  1.  The  mass  ROI  from  an  individual  image  is  used  as  input.  The  mass  is  segmented  from  its  breast  tissue 
background.  Three  types  of  image  features,  namely,  spiculation,  texture,  and  morphological  features  are  extracted  from 
the  segmented  mass  to  describe  the  size,  shape,  margin  characteristics,  and  texture  in  the  breast  tissue  surrounding  the 
mass.  For  mass  classification  using  multiple  images,  the  corresponding  features  extracted  from  the  different  images  are 
averaged.  To  reduce  the  dimensionality  of  the  feature  space,  stepwise  feature  selection  is  performed  to  identify  the  most 
effective  features.  A  feature  classifier  is  then  designed  for  the  differentiation  of  malignant  and  benign  masses  in  the 
multi-dimensional  feature  space. 

For  an  input  ROI  image  from  a  DBT  slice  or  a  PV  image,  the  mass  is  segmented  by  an  active  contour  model  initialized 
with  adaptive  k-means  clustering  followed  by  morphological  filtering.  For  extraction  of  spiculation  features,  a 
spiculation  likelihood  map  is  generated  by  analyzing  the  gradient  directions  around  the  mass  margin  and  spiculation 
measures  are  extracted  from  the  map.10  For  extraction  of  texture  features  in  tissues  around  the  mass,  the  rubber  band 
straightening  transform  (RBST)  is  applied  to  a  band  of  pixels  around  the  segmented  mass  boundary  such  that  the  mass 
boundary  was  in  the  horizontal  direction  and  the  spiculations  are  approximately  in  the  vertical  direction. 1 1  The  gradient 
of  the  RBST  image  is  enhanced  by  Sobel  filtering  in  both  the  horizontal  and  vertical  directions.  The  run-length  statistics 
(RLS)  texture  features  are  then  extracted  from  the  gradient  images.  Morphological  features  including  the  mass  size,  the 
Fourier  descriptor,  and  those  from  the  normalized  radial  length  (NRL)  are  extracted  to  describe  the  mass  shape  and  size. 

For  analysis  using  the  reconstructed  DBT  slices,  we  compared  a  feature  space  obtained  from  features  of  the  central  slice 
alone  and  a  second  feature  space  obtained  by  averaging  the  corresponding  features  from  5  slices  centered  at  the  central 
slice.  For  classification  of  the  masses  using  the  PVs,  we  evaluated  a  feature  space  obtained  from  the  central  PV  alone 
and  another  by  averaging  the  corresponding  features  from  all  1 1  PVs. 

2.3  Feature  Classification 

In  each  feature  space,  a  linear  discriminant  analysis  (LDA)  classifier  with  stepwise  feature  selection  was  trained  and 
tested  using  a  two-loop  leave-one-case-out  resampling  procedure.  In  the  “outer”  leave-one-case-out  loop,  one  case 
including  all  mass  ROIs  from  the  same  patient  was  left  out  in  each  cycle.  Simplex  optimization  was  used  to 
automatically  select  the  best  set  of  parameters  for  a  stepwise  feature  selection  procedure  within  the  training  set  of  (N-l) 
cases  in  that  cycle.  During  the  simplex  optimization  process,  an  “inner”  leave -one-case-out  loop  was  used  to  obtain  an 
area  (Az)  under  the  receiver  operating  characteristic  (ROC)  curve  from  the  (N-l)  validation  cases  as  a  figure-of-merit  for 
guiding  the  search.  The  best  set  of  parameters  was  then  used  to  select  features  from  the  entire  training  set  of  (N-l)  cases 
and  an  LDA  classifier  was  formulated.  The  LDA  was  applied  to  the  left-out  case  to  obtain  a  test  discriminant  score. 
The  same  process  was  repeated  for  all  N  cases  in  the  data  set  in  the  outer  leave-one-case-out  loop.  After  the  two-loop 
resampling  procedure  was  completed,  the  test  scores  for  all  N  cases  were  collected  and  evaluated  by  ROC  analysis.  The 
classification  accuracy  of  the  classifier  was  estimated  as  the  test  Az. 
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The  same  feature  selection  and  classification  procedure  was  performed  in  each  of  the  feature  spaces  as  shown  in  Fig.  1. 
The  test  ROC  curves  and  the  Az  values  were  then  compared  for  the  four  conditions. 


Fig.  1.  Mass  classification  scheme  for  both  approaches  using  the  PVs  as  input  and 
using  the  DBT  slices  as  input. 


3.  RESULTS 

An  example  of  an  ROl  in  the  PV  and  the  DBT  of  a  breast  with  a  spiculated  mass  is  shown  in  Fig.  2.  The  ROIs  from  the 
PV  at  0  degree  and  from  a  SART  reconstructed  slice  through  approximately  the  center  of  the  mass  are  shown. 


(a)  (b) 


Fig.  2.  An  example  of  a  spiculated  mass  in  the  data  set.  (a)  the  ROI  image  in  the  central  projection 
view  (0  degree),  and  (b)  the  ROI  from  the  SART-reconstructed  slice  intersecting  approximately 
the  center  of  the  same  mass.  The  spicules  radiating  from  the  mass  can  be  seen  much  more 
clearly  in  the  DBT  slice. 
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The  classification  of  masses  using  the  reconstructed  DBT  slices  achieved  a  test  Az  of  0.84+0.04  and  0.91+0.03, 
respectively,  in  the  central  DBT  slice  feature  space  and  in  the  5-DBT-slice  feature  space.  The  difference  in  the  Az 
between  the  central  DBT  slice  and  5-DBT-slice  feature  spaces  was  statistically  significant  (p=0.003).  For  classification 
of  masses  using  the  PVs,  the  test  Az  was  0.79+0.04  and  0.83+0.04,  respectively,  in  the  central  PV  feature  space  and  in 
the  11-PV  feature  space.  The  difference  in  the  Az  between  the  central  PV  and  11 -PV  feature  space  did  not  achieve 
statistical  significance  (p — 0. 1 8).  Finally,  the  classification  accuracy  in  the  5-  DBT-slice  feature  space  was  significantly 
better  (p — 0.006)  than  that  in  the  11-PV  feature  space.  The  Az  values  are  summarized  in  Table  1  and  the  ROC  curves 
obtained  under  the  four  conditions  are  compared  in  Fig.  3. 

Table  1.  The  classification  performance  of  the  LDA  classifiers  in  the  four  feature  spaces. 


Az  (PV) 

Az  (DBT) 

Single  PV:  0.79+0.04 

Central  slice:  0.84+0.04 

11  PVs:  0.83+0.04 

5  slices:  0.91+0.03 

Fig.  3.  ROC  curves  of  the  LDA  classifiers  in  the  four  feature  spaces:  (1)  features  extracted 
from  a  single  DBT  slice  intersecting  approximately  the  center  of  the  mass,  (2) 
features  extracted  from  5  DBT  slices  centered  at  the  central  slice  and  averaged,  (3) 
features  extracted  from  a  single  PV  image  at  the  center  of  the  tomosynthesis  scan 
(0  degree  projection  angle),  and  (4)  features  extracted  from  1 1  PVs  centered  at  the 
0  degree  PV  and  averaged. 
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4.  DISCUSSION  AND  CONCLUSIONS 


DBT  is  a  new  imaging  modality  that  holds  the  promise  of  improving  breast  cancer  diagnosis  by  reducing  the  overlapping 
dense  tissue  that  may  obscure  the  characteristics  of  mass  margins.  Development  of  CAD  systems  for  DBT  is  still  at  an 
early  stage.  In  this  preliminary  study,  we  compared  the  effectiveness  of  the  mass  characteristics  obtained  by 
computerized  feature  analysis  on  the  reconstructed  DBT  slices  and  on  the  PV  images  to  differentiate  malignant  and 
benign  masses.  To  our  knowledge,  this  is  the  first  study  that  compared  mass  classification  in  DBT  using  the  two 
different  approaches. 

Our  results  demonstrate  that  the  features  of  breast  lesions  extracted  from  DBT  slices  could  provide  significantly  higher 
accuracy  than  those  from  PV  images  for  classification  of  malignant  and  benign  masses.  Furthermore,  the  feature  spaces 
obtained  by  averaging  corresponding  features  from  several  DBT  slices  or  from  the  PVs  could  provide  significantly 
higher  accuracy  than  those  from  the  central  slice  or  from  the  central  PV  alone,  respectively.  Although  all  image 
information  is  already  contained  in  the  set  of  PV  images,  tomosynthesis  reconstruction  can  combine  the  information 
from  the  multiple  PVs  more  efficiently  than  simply  averaging  the  features  extracted  from  the  individual  PVs.  This 
results  in  the  higher  classification  accuracy  achieved  using  features  from  the  DBT  slices  than  using  features  from  the 
PVs.  For  a  given  approach,  the  average  features  from  multiple  PVs  or  DBT  slices  provided  a  higher  accuracy  than 
features  from  a  single  image  probably  because  of  noise  reduction.  Study  is  underway  to  further  improve  both  approaches 
and  to  evaluate  the  effects  of  combining  the  feature  spaces.  A  larger  data  set  is  also  being  collected  to  improve  the 
training  of  the  CADx  systems. 
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Abstract.  Digital  Breast  Tomosynthesis  (DBT)  Mammography  is  an  emerging 
technique  that  has  the  potential  to  improve  breast  cancer  detection.  In  DBT, 
low-dose  mammograms  are  acquired  at  a  number  of  projection  angles  over  a 
limited  range  and  the  3D  breast  volume  is  reconstructed.  In  this  study,  we 
investigated  the  effect  of  different  distributions  of  projection- view  (PV)  images 
that  included  different  angular  range  and  angular  spacing  on  the  reconstruction 
image  quality.  A  GE  prototype  DBT  system  was  used  to  acquire  a  total  of  21 
PVs  in  3°  increments  over  a  ±30°  range,  from  which  multiple  subsets  containing 
the  same  number  of  1 1  PV  images  were  selected.  A  custom-built  breast 
phantom  and  a  selected  patient  case  were  used  to  evaluate  the  image  quality. 
For  breast  phantom  study,  the  contrast- to- noise  ratio  (CNR),  the  normalized 
line  profiles  of  test  objects,  and  an  artifact  spread  function  (ASF)  were  used  as 
performance  measures  to  compare  the  results  for  the  different  subsets  and  for 
the  full  set.  The  simultaneous  algebraic  reconstruction  technique  (SART)  was 
used  to  reconstruct  the  DBT  under  all  conditions.  Our  results  demonstrated  that 
large  DBT  angular  range  gave  superior  CNR  and  ASF  for  masses  and  less 
interplane  blurring  for  high-density  objects.  Narrow  angular  range  favors  in¬ 
plane  edge  sharpness  for  high-density  objects. 

Keywords:  Digital  Breast  Tomosynthesis  (DBT)  Mammography,  distributions 
of  projection-view  (PV)  images,  simultaneous  algebraic  reconstruction 
technique  (SART). 


1  Introduction 

Digital  breast  tomosynthesis  (DBT)  mammography  is  an  emerging  technique  that  can 
provide  volumetric  information  for  the  breast.  The  breast  volume  is  reconstructed 
from  multiple  low-dose  projection-view  (PV)  images  acquired  at  a  number  of  angles 
in  a  limited  angular  range  [1,2].  The  total  radiation  dose  is  set  to  be  comparable  to 
that  used  in  regular  mammography.  The  depth  resolution  obtained  in  DBT  can  reduce 
the  camouflaging  effect  of  the  overlapping  fibroglandular  breast  tissue  in 
conventional  mammography  and  improve  the  conspicuity  of  subtle  lesions,  and  thus 
potentially  improve  breast  cancer  detection  at  its  early  stage  [3,4]. 


Different  from  computerized  tomography,  DBT  is  a  limited-angle  imaging 
modality.  The  system  design,  including  but  not  limited  to  the  parameters  such  as  the 
number  of  PVs,  the  total  angular  range  and  the  angular  spacing,  needs  systematic 
investigation.  The  reconstruction  image  quality  depends  strongly  on  the  selection  of 
these  and  other  parameters  and  the  reconstruction  algorithm.  In  this  work,  our  purpose 
is  to  investigate  the  effect  of  different  distributions  of  PV  images  on  the  reconstructed 
image  quality  based  on  a  fixed  total  dose.  The  results  will  provide  insight  into  the 
design  of  DBT  systems  and  information  for  optimizing  imaging  conditions. 


2  Methods 

2.1  DBT  System 

A  GE  GEN2  DBT  Mammography  system  was  used  to  acquire  a  total  of  21  PVs  in  3° 
increments  over  a  ±30°  range  [5].  The  imaging  geometry  of  this  DBT  system  is 
schematically  shown  in  Figure  1 .  The  distance  from  the  x-ray  focal  spot  to  the  center 
of  rotation  is  64  cm  and  the  focal-spot-to-detector  distance  is  66  cm.  The  system  has  a 
Csl  phosphor/a:Si  active  matrix  flat  panel  digital  detector  with  a  pixel  size  of  0.1 
mmxO.l  mm.  The  image  acquisition  process  takes  less  than  8  seconds  and  the  digital 
detector  (X-Y  plane)  is  stationary  during  image  acquisition.  The  system  is  operated  at 
an  exposure  level  that  acquires  a  DBT  mammography  at  about  1.5  times  of  the  dose 
of  a  conventional  screen-film  mammogram. 

In  a  preliminary  study  [6],  we  demonstrated  that  the  simultaneous  algebraic 
reconstruction  technique  (SART)  [7]  can  provide  high-quality  tomosynthesized  slices 
while  being  more  efficient  than  statistical  reconstruction  methods  for  specific  system 
configurations  and  imaging  conditions.  In  this  study,  we  employed  SART  for  DBT 
reconstruction.  The  imaged  volume  was  subdivided  into  an  array  of  voxels,  of  which 
the  X  and  Y  dimensions  were  set  to  be  0. 1  mm,  the  same  as  the  detector  pixel  size, 
and  the  Z  dimension  (slice  interval  or  spacing)  1.0  mm.  In  SART,  an  update  of  the  3D 
volume  of  the  attenuation  properties  of  breast  tissue  was  performed  using  individual 
PV  sequentially. 


2.2  Different  PV  Distributions 

Among  many  factors,  the  distribution  of  PV  angles  plays  an  important  role  in  DBT 
system  design.  The  total  angular  range  and  the  angular  sampling  interval  directly 
affect  the  reconstruction  image  quality.  For  low-contrast  objects  such  as  masses,  and 
high-contrast  objects  such  as  calcifications,  the  same  PV  distribution  may  have 
different  influences  on  the  detection  and  characterization. 

To  investigate  the  effect  of  the  distribution  of  PV  angles  on  DBT  image  quality,  six 
subsets,  each  containing  11  PVs,  were  selected  from  the  original  21  PVs.  These  six 
subsets  fell  into  three  categories,  namely,  uniform  group,  non-uniform  central  group, 
and  non-uniform  extreme  group.  The  uniform  group  included  two  conditions, 
representing  a  narrow  and  a  wide  angular  range,  in  which  each  PV  was  spaced  in 
equal  angular  intervals.  The  two  non-uniform  groups  also  contained  two  conditions  in 


each:  the  central  group  had  more  PVs  in  the  middle  region  while  the  extreme  group 
had  more  PVs  in  the  extreme  side  regions.  The  two  conditions  in  each  non-uniform 
group  slightly  differed  in  distribution.  The  total  dose  for  each  subset  was  equal 
because  they  contained  the  same  number  of  PVs.  In  tomosynthesis  reconstruction,  it 
is  known  that  PV  images  acquired  close  to  the  middle  region  are  important  for  the 
spatial  resolution  on  the  X-Y  plane,  while  those  acquired  close  to  the  extreme  side 
region  are  important  for  the  depth  resolution.  The  PV  distributions  in  this  study  were 
chosen  to  evaluate  the  trade-off  among  those  non-uniform  distribution  groups  and  the 
uniform  groups. 

The  x-ray  source  positions  in  terms  of  projection  angles  for  the  six  subsets  are 
listed  in  Table  1  and  schematically  shown  in  Figure  2.  The  six  subsets  are  referred  to 
as  Uniform-Narrow  (UN),  Uniform-Wide  (UW),  Nonuniform-Central-Dense  (NCD), 
Nonuniform-Central-Sparse  (NCS),  Nonuniform-Extreme-Dense  (NED),  and 
Nonuniform-Extreme-Sparse  (NES)  sampling,  respectively.  The  DBT  reconstruction 
using  all  21  PVs  available  is  referred  to  as  the  Full-Set  (FS).  Two  SART  iterations 
were  used  in  the  reconstruction  of  all  subset  conditions,  providing  a  total  of  22 
updates  (the  number  of  updates  in  one  SART  iteration  is  equal  to  the  number  of  PVs 
used  for  reconstruction),  which  was  approximately  equivalent  to  that  of  one  iteration 
when  the  full  set  was  used.  The  relaxation  parameters  were  set  to  be  0.5  and  0.3  for 
the  two  iterations  based  on  a  visual  comparison  of  the  reconstructed  images  of  a 
patient  case  at  different  relaxation  parameter  combinations. 


2.3  Breast  Phantom  and  Figures  of  Merit 

To  quantitatively  evaluate  the  reconstructed  image  quality,  a  custom-built  breast 
phantom  was  used  which  contained  test  objects  on  a  thin  feature  layer  sandwiched 
between  five  1-cm-thick  homogeneous  Lucite  plates  below  and  another  1-cm-thick 
plate  above.  The  test  objects  included  layers  of  thin  circular  aluminum  foils  to 
simulate  masses  with  different  contrasts  and  high-contrast  metal  wires  to  simulate 
metal  markers.  The  wire  was  placed  approximately  perpendicular  to  the  chest  wall 
edge  in  order  to  evaluate  the  DBT  blurring  along  the  x-ray  source  moving  direction 
which  was  parallel  to  the  chest-wall  edge.  The  contrast-to-noise  ratio  (CNR),  the 
normalized  line  profiles  of  test  objects,  and  an  artifact  spread  function  (ASF)  were 
used  as  performance  measures  to  compare  the  results  of  different  subset  conditions 
and  that  of  using  the  full  set.  The  CNR  of  a  simulated  mass  was  defined  as  the 
difference  in  the  average  pixel  intensities  between  the  mass  and  its  local  image 
background,  divided  by  the  standard  deviation  of  pixel  intensity  in  the  same 
background.  The  ASF  of  the  mass  was  defined  as  the  ratio  of  the  CNR  values 
between  a  given  off-focus  layer  and  the  feature  layer,  and  the  ROIs  extracted  for  the 
calculation  of  CNR  values  on  different  layers  were  located  in  the  same  positions  in 
terms  of  the  X  and  Y  coordinates.  Thus,  the  ASF  quantified  the  inter-plane  artifacts 
along  the  Z-direction.  The  line  profile  of  a  mass  was  extracted  from  the  pixel 
intensities  along  the  scan  direction  of  the  x-ray  source  (Y-axis)  with  the  removal  of 
the  mean  pixel  intensity.  Similarly,  to  facilitate  the  comparison  of  relative  sharpness, 
line  profiles  of  the  wires  were  extracted  along  both  the  Y  and  Z  axes,  and  the  resulting 
maximum  intensity  was  normalized  to  1  after  removal  of  the  mean  pixel  intensity. 


3  Results 


The  CNR  values  of  selected  simulated  masses  with  different  contrasts  were  listed  in 
Table  2  for  all  imaging  conditions.  The  CNR  values  are  sorted  approximately  in  an 
ascending  order  with  the  corresponding  total  angular  range,  and  the  number  of  PVs 
outside  ±6°  is  also  listed  in  order  to  illustrate  the  sampling  density  in  the  extreme  side 
regions.  The  results  suggested  that  the  mass  CNR  value  depended  primarily  on  the 
total  angular  range  and  to  some  extent  on  the  number  of  PVs  acquired  at  large 
projection  angles.  Among  the  six  subsets,  the  UW  and  NED  sampling  that  have  more 
PVs  in  the  extreme  side  regions  and  the  largest  angular  range  provided  the  best  CNRs, 
probably  because  the  large-angle  PVs  can  better  separate  overlapping  in  the  depth 
direction  and  thus  enhance  the  object  contrast  in  the  reconstruction.  All  CNR  results 
from  the  subsets  were  inferior  to  those  from  the  full  set. 

The  line  profiles  and  the  ASF  curves  of  the  low-contrast  mass  were  shown  in 
Figures  3  and  4,  respectively.  The  results  of  the  high-contrast  mass  were  not  shown 
but  have  similar  behavior.  There  is  no  obvious  difference  among  the  normalized  line 
profile  of  masses  from  the  six  subsets,  probably  because  the  line  profiles  of  relatively 
low  contrast  objects  are  not  sensitive  to  in-plane  blurring  due  to  the  subset  sampling. 
The  UN  sampling,  which  had  the  narrowest  angular  range  (±15°)  among  all 
conditions,  gave  much  worse  ASF  behavior  than  others.  Among  the  rest,  the  NCD 
sampling,  which  had  slightly  narrower  angular  range  (±27°)  than  others  of  full  range 
(±30°),  gave  slightly  worse  ASF.  In  other  words,  interplane  blurring  of  mass 
depended  strongly  on  the  total  angular  range.  The  NED  sampling,  which  had  all  10 
PVs  evenly  distributed  on  extreme  sides  in  addition  to  one  in  the  middle,  yielded  the 
best  ASF  or  least  interplane  blurring,  even  slightly  better  than  the  full  set  result.  This 
is  probably  because  the  PVs  in-between,  especially  those  close  to  the  middle  region, 
contribute  more  blurring  to  the  out-of-focus  planes  than  the  extreme  PVs. 

The  normalized  in-plane  and  interplane  line  profiles  of  the  high-contrast  wire  were 
shown  in  Figures  5  and  6,  respectively.  Different  from  the  simulated  low-contrast 
mass,  the  UN  result  provided  the  best  in-plane  sharpness  while  the  NED  sampling 
worst.  This  is  because  the  extreme  side  PV  distributions  did  not  provide  enough 
information  to  reconstruct  the  high-contrast  wire  on  the  in-focus  plane,  thus  reducing 
the  sharpness.  On  the  other  hand,  based  on  the  similar  reason  we  discussed  about 
simulated  masses,  the  NED  sampling  provided  the  narrowest  ASF  and  thus  least  inter¬ 
plane  blurring  for  the  high-contrast  wire. 

In  addition  to  the  phantom  study,  a  DBT  patient  case  was  reconstructed  with  the 
six  subset  sampling  and  the  full-set  of  PVs.  An  ROI  containing  a  spiculated  mass  was 
extracted  for  comparison.  Figure  7  shows  the  ROIs  from  the  middle  PV  image  and 
from  a  DBT  slice  intersecting  approximately  the  center  of  the  mass  reconstructed 
under  each  of  the  conditions.  It  can  be  seen  that  on  the  middle  PV  image,  the  mass 
boundary  and  the  spiculations  appear  to  be  blurred  due  to  the  overlapping  tissues.  The 
full-set  reconstruction  results  in  the  best  image  quality  with  sharp  speculations,  and 
the  breast  parenchyma  is  reasonably  free  of  artifacts.  In  comparison,  the  UN  subset 
yields  relatively  sharp  spiculations  but  there  are  additional  structures  on  the  image 
that  may  be  caused  by  overlapping  tissues  from  the  reduced  depth  resolution.  The 
spiculations  and  edges  from  the  UW  subset  are  slightly  more  blurred  than  those  from 
the  UN  subset  but  the  background  exhibits  less  structures.  The  NCD,  NCS  and  NES 


subsets  provide  comparable  spiculations.  The  NED  subset  results  in  the  worst 
spiculation  details  and  most  blurred  mass  margins.  The  textural  artifacts  in  the 
background  tissue  from  the  NCD.  NCE,  NED,  NES  subsets  are  worse  than  those  from 
the  UN  and  UW  subsets  in  some  regions. 


4  Discussion 

DBT  mammography  can  provide  quasi- 3D  information  of  the  breast  in  contrast  to 
conventional  mammography.  With  a  breast  phantom  study,  we  compared  results  of 
six  subset  samplings  of  1 1  PVs  out  of  the  21  original  PVs  and  the  results  of  the  full 
set  of  21  PVs.  Preliminary  results  demonstrated  that  the  larger  the  DBT  angular  range, 
the  better  the  CNR  and  ASF  of  simulated  masses.  In  case  of  equivalent  angular  range, 
PVs  acquired  at  large  projection  angles  had  stronger  effect  on  depth  resolution  than 
those  close  to  the  middle  region.  For  a  high-contrast  wire,  a  narrower  angular  range 
led  to  shaper  in-plane  edge  and  the  PVs  acquired  at  large  projection  angles  reduced 
the  blurring  of  inter-plane  profile.  No  subset  was  superior  to  the  others  in  all 
performance  measures,  indicating  the  inevitable  tradeoff  between  in-plane  resolution 
and  depth  resolution.  The  tradeoff  is  less  clear  for  the  DBT  images  of  a  spiculated 
breast  mass.  The  overlapping  tissue  causes  structural  textures  on  the  reconstructed 
slices.  The  full  set  reconstruction  demonstrates  a  clear  advantage  in  terms  of  in-plane 
and  depth  resolutions.  The  superior  image  quality  can  be  partly  attributed  to  the 
higher  dose  with  the  use  of  all  21  PVs.  The  comparison  of  the  PV  distribution  in  the 
full  set  with  those  of  the  other  six  subsets  is  therefore  not  completely  fair  unless  we 
can  keep  the  total  dose  of  the  set  of  PVs  at  a  constant  level.  Due  to  restrictions  on  the 
unit  required  for  other  clinical  evaluations,  our  DBT  system  is  not  enabled  to  acquire 
PV  images  of  different  dose  distributions.  However,  it  has  been  suggested  that  higher 
dose  acquisition  in  some  specific  PVs  than  others  might  improve  image  quality.  We 
will  include  this  non-uniform  dose  approach  in  future  studies  if  advanced  system 
operating  modes  become  available  or  by  using  computer  simulations. 
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Table  1.  Projection  angles  (in  degrees)  of  the  six  different  sampling  of  PVs. 


Uniform-Narrow  (UN) 

±15 

±12 

±9 

±6 

±3 

0 

Uniform-Wide  (UW) 

±30 

±24 

±18 

±12 

±6 

0 

Nonuniform-Central-Dense  (NCD) 

±27 

±18 

±9 

±6 

±3 

0 

Nonuniform-Central-Sparse  (NCS) 

±30 

±21 

±12 

±6 

±3 

0 

Nonuniform-Extreme-Sparse  (NES) 

±30 

±27 

±24 

±6 

±3 

0 

Nonuniform-Extreme-Dense  (NED) 

±30 

±27 

±24 

±21 

±18 

0 

Table  2.  The  contrast-to-noise 

ratio  (CNR)  results  of  the  selected  masses. 

UN 

NCD 

NES 

NCS 

UW 

NED 

FS 

Angular  range  +15 

±27 

±30 

#  PVs  outside  ±6°  6 

6 

6 

6 

8 

10 

16 

Low-contrast  mass  2.85 

3.56 

3.68 

3.88 

4.13 

4.05 

4.74 

High-contrast  mass  6.36 

7.37 

7.51 

7.85 

8.27 

8.28 

9.47 

Fig.  1.  Geometry  of  the  GE  GEN2 
digital  breast  tomosynthesis 
mammography  system. 
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Fig.  2.  Six  different  distributions  of  PV  images. 


Fig.  3.  Normalized  in-plane  line  profiles  of 
the  selected  low-contrast  mass  for  different 
subsets  of  PVs  and  the  full  set. 


Distance  from  feature  layer  (mm) 

Fig.  4.  ASF  of  the  selected  low-contrast 
mass  with  different  subsets  and  the  full  set. 
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Fig.  5.  Normalized  in-plane  line  profiles  of  Fig.  6.  Normalized  inter-plane  line  profiles 
the  selected  wire.  along  Z-axis  of  the  selected  wire. 


UW  NCS  NED 


Fig.  7.  Spiculated  mass  of  the  selected  patient  DBT  case  on  projection  view  image  (0  degree) 
and  reconstructions  under  the  full  set  and  the  six  subset  sampling  conditions. 


